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SECTION 1 


INTRODUCTION 


This report is meant to supplement Draper Laboratory Report R-1113, 

111 * - 

Optimal Solar Sail Planetocentric Trajectories . it contains informa- 

tion about the computer code, called SUNSPOT (for SUN-Sail Program for 
Optimal Trajectories), which was used to generate the results in the 
report, and should allow one to use this code to produce additional 
solar sail trajectories. The principal contents are a description of 
the input and output and brief descriptions of the subroutines and list- 
ings of all the subroutines. The code is functionally similar to 
SECKSPOT, a solar electric orbit transfer program. (2,3) 

The program calculates the solution to a two-point boundary- value 
problem which arises from the application of optimal control theory. 

A general flow chart is shown in Figure 1. Figure 2 is a system diagram 
which illustrates the connection of various subroutines which make up - 
the deck, The subroutines above call those below. 

Total size of the object deck compiled on an Amdahl 470 in Fortran 

H is approximately 59,000 bytes. The size of the individual subroutines 

is shown in Table 1. A Fortran Scientific Subroutine Package (SSP) rou- 

(4) 

tine called DRTNI is also needed, but not shown in the table. 

A plotting capability was developed at Draper Laboratory. Data 
can be written by the output subroutine and stored for later plotting. 
Since plotting capabilities vary from facility to facility, that func- 
tion is not discussed in this report. 



Superscript numerals refer to similarly numbered items in the 
List of References. 
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. m able i. Program subroutine object deck size 
(thousands of bytes! . 


CONE 

0.9 

OUTP 

4.8 

CONTL 

0.6 

OUTPC 

2.2 

DCUBIC 

1.2 

PLANET 

1.5 

DCROUT 

2.5 

PRTN 

0.6 

DQRTIC 

1.3 

QUAD 4 

1.0 

EVALMP 

5.2 

QUAD 8 

1.2 

FCT 

3.4 

QUAD 1 6 

1.8 

FUNCT 

2.5 

QUAD 32 

2.9 

INPUT 

9.5 

RK4 

1.7 

ITER (NR) 

3.2 

SHADOW 

3.5 

ITER (MODNR) 

3.9 

SUN 

1 .4 

OBLATE 

2.0 

TRAJ 

2.1 
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SECTION 2 
INPUT 

The quantities discussed in this section are all read by the sub- 
program INPUT. Unless otherwise indicated each value is read on a 
separate line, real variables in fixed format (F25, 15) , integer variables 
beginning with i, j, k, 1, m, n are read in format, 12. The inputs for 
the deck are listed in this section along with a brief description and 
nominal values, if any. 

IPNUM planet flag: 1 = Mercury, 2 = Venus, 3 = Earth, 4 = Mars 

(format, II) 

Initial orbit 

W(l) (km) semimajor axis 

W(2) eccentricity 

W(3) (degrees) inclination 

W{4) (degrees) longitude of ascending node 

W(5) (degrees) argument of pericenter 

Initial guesses 

ZL0(1) A, adjoint to semimajor axis 

ZLQ(2) A^, adjoint to orbital element h 

ZL0(3) , adjoint to orbital element k 

ZLO (4) A , adjoint to orbital element p 

r 

ZLO (5) Ag, adjoint to orbital element q 

The desired final orbit 

WF(1) (km) semimajor axis 

WF(2) eccentricity, not used if NIP = 3 

RF(3) (deg) inclination, not used if NIP = 3 

5 



n- 





(deg) longitude of ascending node, i;oc used if NOP - 2 or 3 
(deg) argument of pericenter, not used if NOP = 2 or 3 

(days) guess for final time 
2 

(rran/s~) characteristic acceleration. 

Julian date at initial time 

The following input nay be read or, optionally, left at nominal 
values, IRDFLG is read followed by the addition input or operations and 
then IRDFLG is read again until IRDFIjG = 01 and input is ended. 


WF ( 4 ) 
WF (5) 
TF2 
AC 
TL 


IRDFLG NOMINAL 


1 End of input 


2 

IPR 

print flag (for intermediate 
traj ectories) 

0 

3 

NIMAX 

maximum number of iterations (if 0, 
bypass iteration to print time 
history) 

20 

4 

TFMAX2 

(days) maximum flight time 

500. 

5 

DT2 

(days) time step for integrator 

10. 

6 

UEB 

upper error bound for integrator 
(not used with RK4 ) 

1.D10 

7 

EV1 

error weights for integrator 
(5D6.1, 2 lines) (not used with RK4) 

1., 1,1, 1,1,0 . 

8 

UTKM 

equatorial radius (km) 

Set in PLANET 

9 

GM 

(km3/s2) gravitational constant 

Set in PLANET 

10 

NO? 

= 1, five orbital elements specified 
at TF 

1 



= 2, three orbital elements specified 
at TF 




= 3, semimajor axis specified at TF 


11 

Sets 

oblateness 

AJ2, = 0 .DO 

Set in PLANET 

12 

STEP 

step size for numerical differentia- 
tion in ITER, 6 dim., sixth element 
for time variation of Hamiltonian 

l.D-2 


KSTEP 

= 0, step as fraction in ITER 

1 


= 1, step as constant in ITER 
(except STEP (6)) 
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0 


> 


ISON = 0, shadow effect off 

= 1, shadow effect on 
14 ISUN = 0, sun distance effect on thrust off 0 


= 1, effect on 


15 

CD (1) 

coefficients in 

0.5 


CD ( 2 ) 

acceleration factor 

0.5 


CD (3) 


0.0 

16 

IE 

= 1, if equatorial frame, 

1 


solar motion 

= 2, if equatorial frame, 
no solar motion 

= 3, if e c 1 i p t i c - frame, solar notion 

= 4, if ecliptic frame, no solar motion 

= 5, if planetary frame, solar motion 

= 6, if planetary frame, no solar 
motion 
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FLIM 

norm limit in iteration routine 
(format, D6.1) 

l.D—6 

18 

XNU 

pericenter penalty function 
weighting 'factor 

0.0 

19 

NORB 

= 0, no orbit print 

= 1,..., 999, orbit print on 

NORB points of an orbit (format, 13) 

0 

20 

Sets flag 

so that plot data is stored 


21 

IQ 

orbit divided into IQseparate 
quadrature intervals (1-10) 

2 


. J 
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OUTPUT 


Most of the output is self-explanatory and a look at an example 
will familiarize the user with it. There are certain basic groups of 
output. The first is the printing of the read-in initial data and a 
few internally set constants. Normally this will be followed by output 
from the iterator. After convergence, a summary of characteristics of 
the converged trajectory is printed. Finally, a tine history of the 
converged trajectory will be printed. Usually, even if convergence 
was unsuccessful, a time history of the last trajectory to be calculated 
will be printed. 

The printing of the initial data should be understandable. There 
are a few abbreviations used. 


A 

semimajor 

axis 



E 

eccentricity 



I 

inclination 



I£N ASC NODE 

longitude 

of ascending node 

ARG PERIC 

argument 

of pericenter 

MS 

meters/second 



P.R. 

planet radii 



UTKM 

internal 

units 

t o 

ki lomete.rs 

UTS 

internal 

units 

t o 

seconds 

UTD 

internal 

units 

t o 

days 

UTKG 

internal 

units 

to kilograms 

UTKW 

internal 

units 

t o 

kilowatts 

UTMS2 

internal 

units 

t o 

meters/second2 
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After the initial input print, the iteration begins. The iteration 
number (ITER NO.) and the total number of calls to TRAJ are printed fol- 
lowed by X, the iteration parameters (XLO), then Y, the error in the 
final conditions. The final conditions are the final values of a, h, k, 
p, q, H, if NOP = 1. If NOP = 2 the final conditions are a, e, tan 
X , A^, H unless the final eccentricity or inclination is zero. 

If eccentricity is zero the second condition is h and the fourth k; if 
inclination is zero, the third is p and the fifth is q. If NOP = 3, the 

final conditions are a. A, , A, , A. A , H. Then the final time (TF) is 

** K p ■ 

printed in internal units, followed by, FO, the sum of the squares of 
the errors in the final conditions. For convergence this value must be 
less than FLIM, the "norm limit in ITER". In order to calculate the 
partial derivative matrix or sensitivity matrix, the nominal values of 
"X" are changed slightly by inputted amounts; these perturbed values Of 
X(X(I) + DX(1)) are next printed followed by the corresponding Y. Then 
the partial derivative matrix is printed as well as its determinent. 

This matrix is inverted and premultiplies the error vector to obtain the 
changes in the X's, DELX:S, which are next printed. 

A new trajectory is calculated and the sun of the squares of the 
errors in the final conditions is printed (FI). If this is smaller 
than FO, a new iteration begins; if it is larger than FO, the DELX:S 
are halved and printed. This continues until FI < FO or until a certain 
number of halvings. What follows depends on how well the method converges 
and 'on whether the Newton-Raphson or modified Newton-Raphson subprogram 
is used. Further output is basically permutations of the above, termin- 
ating with convergence or a message indicating lack of success. 

After exit from the iteration, a summary of characteristics of the 
last trajectory (the optimal, if convergence was successful) is printed. 
Included are the actual final orbital elements, the error in the final 
orbital elements, the values for the iteration parameters, and the final 
time. 

Next is printed a time history of the final (optimal if convergence 
was successful) trajectory. If NIMAX = 00, then a time history is printed 
immediately following printing of the input data, bypassing the iteration 
routine, and summary print. 

The low thrust trajectory time history at each time step is printed. 
First is printed TIME in various units. The time step number is also 
printed. Next are printed the equinoctial orbital elements (a, h, k, p, q) 
followed by the classical orbital elements (a, e, i, SI, w) . Next is the 
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costate, then the state derivative, then the costate derivative and then 
the value of the Hamiltonian, the period (hours), pericenter and apocenter 
(km). The characteristic acceleration divided by the sun distance squared, 
and the value of C^* If shadowing is included the time spent in shadow 
is printed in hours and as a fraction of the period (if the orbit passed 
through shadow). The entry angle (°) and exit 'angle (°) are printed. 

If NORB was set to a nonzero value then there is additional output. 

The optionally printed output for an orbit at each time step includes 
the sun direction in the equinoctial coordinate frame and the planet-sun 
distance in units of the planets semimajor axis. Several spacecraft 
parameters are printed at a number of points on the orbit (the number 
of points = NORB). The points are at equal units of eccentric longitude. 

At each point the spacecraft location is printed including the eccentric 
longitude and the £ and g_ components X-^ and in units of planet radii. 

The following spacecraft parameters are printed with the spacing indi- 
cated here: 


PRIMER CONE A THRUST CONE A PSIV ACC FACT TAV RATIO 

UF UG UW 

All angles are in degrees. The meaning of these terms follows, Equa- 
tion numbers refer to Reference 1. 


PRIMER CONE A the primer vector cone angle, 3, as i n Eq, (3.50) 
THRUST CONE A the force vector cone angle, 0 , as in Eq. (3.27) 


PSIV 


ACC FACTOR 
TAV RATIO 
UF, UG, UW 


the force vector clock angle, ip, as in Eq. (3.52), 
except that \p has the opposite sign, where the 

reference vector is the velocity vector 

\ 

the value of the acceleration factor, c, + c- cos 20 

+ cos 49 as in Eq. ( 3 . 26 ) ' 

the ratio of the photon force with the local gravita- 
tional force 

A A A 

the f, g, and w components, respectively, of the thrust 
acceleration unit vector, in the equinoctial coordinate 
system 


A number of error messages are scattered through the code. A few 
will be mentioned here. Several, in INPUT, call attention to bad input 
data. For bad input data following an IRDFLG value, a message, 

IRDFLG = (number], is printed. In some cases additional information is 
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given. Wien shadowing is included, a message, ISHAD = 1, indicates that 
only one shadow crossing was found. This usually arises from small 
numerical inaccuracies in solving the quartic equation and can usually 
be ignored unless a number of orbits have intersected the planet's 
surface. 
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SECTION 4 


COMMON AREAS 


There are several common areas which are shared by several of the 
subroutines. Table 2 lists these common areas and indicates which sub- 
routines share them. Those subroutines which have no common areas are 
not included, Following is a list of common areas, the variables in 
these areas, the definition of variables, and a list of subroutines 
which contain the particular common areas. Equation numbers refer to 
Reference 1. 


Common area: 
A 

AMU 

PI 


A/A, AMU, PI 

Characteristic acceleration divided by semimajor 
axis of planet squared 

Gravitational constant, p 

TT 


Shared by; 

FCT, FUNCT, INPUT, OUTP, TRAJ 


Common area: 
CD (1) 
CD (2) 
CD (3) 
CD (4) 


CCOM/CD ( 4 ) 

Cf , C^ coefficients in acceleration 

factor, Eq. (3.3) 

Cutoff value of 8, Eq. (3.7) 


Shared by : 

FCT, INPUT, OUTP, OUTPC, CONE 


Common area: 
CB 
SB 

THETA 

CTHETA 

STHETA 


CON/CB, SB, THETA, CTHETA, STHETA, C2T, C4T, S2T, S4T 
Cosine of primer vector cone angle 
Sine of primer vector cone angle 
Force vector cone angle 
Cosine of force vector cone angle 
Sine of force vector cone angle 
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Table 2. Common areas. 


SUBROUTINES 
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C2T Cosine of two times the force vector cone angle 

c4T Cosine of four times the force vector cone angle 

S2T Sine of two times the force vector cone angle 

S4T Sine of four times the force vector cone angle 

Shared by: 

FCT, OUIP, CONE 

Common area: DY/DYDT(6) 

DYDT The partial of the final conditions with respect to time 

Shared by : 

ITER, TRAJ 

Common area: ELEM/ZPO (5) , ZPF(5) 

ZPO The initial state 

ZPF The desired final equinoctial orbital elements 

Shared by: 

INPUT, OUIPC, TRAJ 

Common area; F/FLIM, KSTEP 

FLIM The requirement on the limit of the norm, of the errors 

intheiterator. 

KSTEP Flag indicates whether STEP refers to a fixed increment 
or a fractional increment 

Shared by : 

INPUT, ITER 

Common area: INT/IPR, IDIM, IDIM2 , NIMAX 

IPR Print flag used in OUTP and OUTHX 

IDIM Equal to the dimension of the state plus costate (10) 

IDIM2 Equal to dimension of the state (5) 

NIMAX The maximum number of iterations of the iterator 

Shared by: 

CONTL, INPUT, ITER, OUTP, TRAJ 
Common area: IS/ISUN, ISON 

ISUN Flag indicating if solar distance effect is included 

ISON Flag indicating if shadowing is included 

Shared Jby : 

FUNCT, INPUT, OUIP 


1 4 
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Common area: JD/TL 

TL Julian data at launch 

Shared by : 

PLANET, INPUT 

Common area: J2/AJ2 

AJ2 The oblateness coefficient 

Shared by : 

FUNCT, INPUT, PLANET 

Common area: NU/XNU 

XNU Weighting constant in penalty function, Eq. (3.73) 

Shared by : 

FUNCT, INPUT, TRAJ 

Common area: ORBIT/NORB 

NORB Indicates if there should be orbit print and the number 

of points printed 

Shared by : 

INPUT, OUTP 

Common area: ORBITl/VEC ( 3) , PA 

VEC The thrust direction vector 

PA Accelerat ion factor, Eq.. (3.3) 

Shared by: 

FCT, OUTP 

Common area: 0RBIT2/X1, Y1 , PR(2, 2), X1D0T, Y1D0T 

Xl ? component of spacecraft location in orbit (Eq. 3.16, 

Reference 1) 

A # _ 

Y1 g component of spacecraft location in orbit (Eq. 3.17, 

Reference 1) 

PR(2, 2) Partials of XI and Y1 with respect to h and k (Table B-2, 
Reference 1) 

X1D0T f component of velocity, Eq. (3.18) 

Y1D0T g component of velocity, Eq. (3.19) 

Shared by : 

EVALMP, FCT, OUTP 
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Common area: PLNTS/IPNUM 

IPNUN Integer indicating planet number 

Shared by: 

INPUT, PLANET 

Common area: PLOT/IPL 

IPL Flag indicates if data should be written and stored 

for later plotting 

Shared by: 

INPUT, OUIP 


Common area: 
IQ 


Q/IQ 

The number of segments into which the orbit is divided, 
so that the quadrature is called for each segment to 
calculate the thrust effect. 


Shared by: 

FUNCT, INPUT 

Common area: SHAD/FEN, FEX, DFEN{5) , DFEX(5) , ISHAD 

FEN The eccentric longitude at shadow entry 

FEX The eccentric longitude at shadow exit 

DFEN{5) The partial of FEN with respect to the equinoctial 
orbital elements 

DFEX(5) The partial of FEX with respect to the equinoctial 
orbital elements 

ISHAD A flag indicating whether or not a particular orbit 
passed through the planet's shadow 

Shared by : 

FUNCT, OUIP, SfflW 

Common area: S0L/RSUN{3) , RS 

RSLN The unit vector pointing toward the sim from the planet 

RS The distance from planet to sun divided by semimajor-axis 

Shared by : 

FCT, FUNCT, OUTP, SHADOW, SUN 


Common area: T/TF, TO 

TF Flighttime 


TO 


The initial time (measured from launch time) 


Shared by: 

INPUT, ITER, OUIPC, PRTN, TRAJ 
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Common axea: TC/NOP 

NOP Flag indicates whether one or three final or five 

final orbital elements are specified 

Shared by : 

INPUT, OUTPC, TRAJ 

Common area: TERRA/ AE , EC, W, ENE, AN, ECLMTX, EQUMTX 

AE Planet's semimajor axis (in A.U.) 

EC Eccentricity of planet's orbit 

W Planet's argument of perihelion 

ENE Planet's mean orbital motion 

AN The mean anomaly of the planet at the time of launch 

ECLMTX Matrix Cor conversion to ecliptic coordinate frame 

EQUMTX Matrix for conversion to equatorial frame 

Shared by: 

PLANET, SUN, INPUT 

Ccrnmcn area: TRA/TFMAX, DT, UEB, EW(10) 

UMAX The maximum time of flight allowed 

DT The time step used by the differential equation integrator 

UEB The upper error bound (not used by RK4) 

BV Error weights (not used by RK4) 

Shared by: 

INPUT, TRAJ 

Common area: UNITS/UTS, UTM, UTH, UTD, UTKM, DTR, UTKG, UTKW, UIMS?. 

UTS Internal, units to seconds 

UIM Internal units to minutes 

UTH Internal units to hours 

UTD Internal units to days 

UTKM Internal units to kilometers 

DIR Degrees to radians 

UIRG Internal units to kilograms 

UTKW Internal units to kilowatts 

2 

UEVE2 . Internal units to meters/second 
Shared by 

PLANET, INPUT, QUIP, OUIPC 
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Common area : 
WF 


WF/WF ( 5 ) 

The desired final classical orbital elements 


Shared by : 

INPUT, OUTPC 


Common area: XMMM/ZLO (5) , STEP(6), ZERF(6) 


ZLO The iteration parameters 

STEP Step size used to numerically evaluate partial 

derivative (or sensitivity) matrix in iterator 

ZERF The error in the final conditions 


Shared by : 

INPUT, ITER, OUTPC, PRTN, TRAJ 


Common area: 
ZF 
DZ 


Z/ZF (10) , DZ (10) 

The state and costate 

The derivative of the state and costate 


Shared by: 

OUTPC, TRAJ 
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SECTION 5 


SUBROUTINE DESCRIPTION AND LISTINGS 


This section contains descriptions of all subroutines including 

input and output, common areas, subroutines which are called or called 

by. A listing follows the description. A variable in the argument list 

or the common areas is underlined once if it is output, twice if input 

or three times if both. Aliases are sometimes given in parentheses 

following the subroutine neme. In addition to the subroutines contained 

(4) 

in this section, an IBM Scientific Subroutine, DRTN1 is needed. 

Equation numbers refer to the final report, Reference 1. 
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Subroutine CONE 


Description: 

Calculates JU1 and its derivative (see Eq. (3.60) and (3.64)) 
30 

Used to calculate cone angle by the Newton method. 


Argument List: 


THETA, 

FT, DERFT 




THETA 

Iterative 

value 

of 

thrust cone angle 

FT 

Iterative 

value 

of 

ft 

DERFT 

Iterative 

value 

of 

a-^H- 

AO 4 - 


Common Areas: 

CCOM/ CD ( 4 ) 

C °N/CB, SB, DUB (7) 


Called by: 
DRTNI 
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UOU U U O O O 


- CCHE 

2 VERSION OF CORE TO BE USED WITH NEWTON ITERATOR 
CALLED E I DRTNI 

CALCULATES DH/DTHETA AND DERIVATIVE 

SUBROUTINE CONE (THETA , FT , DERFT) 

IMPLICIT REAL*8 (A-H.O-Z) 

COHKOH /CCOK/CD (4) 

CORHOH /CON/CB,SE,DUB(7) 


CTHETA= DCOS (THETA) 

STHETA= D S IN (THETA) 

C2T= DCOS (2.D0*THETA) 

C4T= DCOS (4. DO* THETA) 

S2T= DSIN (2«D0*THETA) 

S4T= DSIN (4. D0*THETA) 

C 

B= CD (1) +CD (2) *C2T+CD (3) *C4T 
BP= ~2. D0*CD (2) *S2T-4.D0*CD(3)*S4T 
BPP= -4. D0*CD (2)*C2T- 16. D0*CD (3)*C4T 
D= CTHETA*CB+STHETA*SB 
DP= -STHETA*CB+CT BET A* SB 
DPP= -D 
C 

FT= BP*D4B*DP 

DEHFT= 2.D0*BP*BP+BPP*D+B*DPP 


RETURN 

END 


0000001 0 
00000 020 
00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00000090 
00000100 
00000110 
00000 120 
00000130 
00000140 
00000150 
00000 160 
00000170 
00000180 
00000190 
00000200 
000 00210 
00000220 
00000230 
00000240 
00000250 
00000260 
000 0027 0 
00000280 
00000290 
00000300 
00000310 
00000320 
00000330 
00000340 
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Subroutine CONTL 


Description: 

This is the main "driver” program. 


Common Areas: 

INT/ IPR , IDIM, IDIM2, NIMAX 

Calls S ubroutines : 

INPUT, ITER, OUTPC, TRAJ 



oort noo ooooon 


CCNTL/CONTLSS 
C SOLAR SAIL 

C THIS IS THE MAID CONTROLLING PROGRAM FOR FINDING THE 
C OPTIMAL TRAJECTORY FOE PLANETOCERTRIC SOLAR SAIL 
C SPIRAL, ESCAPE OR CAPTURE. 

C 

c 

c 

c 

IMPLICIT REAL*8 (A-H,0-Z) , INTEGER (I-N) 

COMMON /INT/IPR, IDIK, IDIM2, NIMAX 
EXTERNAL TRAJ, PRTH 


ALL SETTING OF CONSTANTS AND INITIAL READ AND WHITES 
ARE DONE BY THE SUEEOUTINES INPUT, PLANETS 

CBLL INPUT 

IF (NIMAX. EQ-0) GC TO 10 
WRITE: (6,7 001) 

THE ITERATOR ROUTINE SOLVES THE 2PBVP FOR TEE OPTIMAL TRAJECTORY 

CALL ITE R (KOUNT , HI, TRAJ, PSTN) 

IF (NT .EQ. 9999) WRITE (6,1002) 

WRITE (6,1003) 

OUTPC PRIWTS A SUMMARY OF THE OPTIYAL TRAJECTORY VALUES 
CALL OUTPC (KOUNT) 

TIME AISTORY OF THE OPTIMAL TRAJECTORY IS CALCULATED BPD PRINTED 

10 RRITE (6,1004) 

IPR= 1000000 
CALL TRAJ 

IF (HI. EQ. 9999) VSIT2 (6,1002) 

IF (NI. BE. 9999) YRITE (6,1005) 

STOP 

1001 PORKAT (18H1 ITERATION BEGINS) 

1002 FORMAT (29HO OPTIMIZATION NOT SUCCESSFUL) 

1003 FORMAT (43H1 CONVERGED VALUES FOR OPTIMIZED TRAJECTORY) 

1004 PORKAT (36R1 TIKE HISTORY OF OPTIMAL TRAJECTORY) 

1005 FORMAT (30HO PROGRAM HAS RUN SUCCESSFULLY) 

END 


00000010 
00000020 
00000030 
00000040 
00000050 
00000060 
00000070 
00030090 
00003090 
00000100 
00000 110 
00000120 
00000 130 
00000 140 
03000150 
00000160 
00000170 
00000 180 
000001 90 

oaooo 200 
00000 210 
00000220 
00000230 
00000240 
00000250 
OOOOO 260 
00000270 
00000 280 
00000290 
00000300 
03000310 
00000320 
00000330 
OOOOO 340 
00000350 
00000360 
00000370 
00000380 
00000390 
00000400 
000004 1 0 
00000420 
00000430 
00000440 
00000450 
00000460 
00000470 
00000480 
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Subroutine DCROUT (DCROUTIO) 


Description: 

This has a square matrix and a vector as input. It inverts the 
matrix and premultiplies the vector by the inverse. 


Argument List: 

AA, g, D, E E§, gj, MD 


AA The matrix which is to be inverted 

R As input, the vector to be multiplied by the 

inverted matrix, as output, the resultant vector 

t> The determinant of AA 

E? S An input, a small quantity which is used to check 

for a singular matrix 

NI Dimension of the matrix and vector 

M Flag set to I as input 

IND Flag, if 0, matrix is singular 


Called by : 


ITER 



BCR 0 UT/D C BOUT 1 0 6 DTH. 

SUEHODTINE DCROUT {l&, R,D,EPS,NI,M,IND) 
DOUBLE PRECISION A, R , D , EPS, T, S, P, DT , ft A 
DIRENSIOH A (6,6) ,R (6,1) ,A A (6,6) 

5 IND=0 
N=NI 

DO 6 1=1, N 
DO 6 J=1,N 

6 A (I,J)=AA(I,J) 

IF(«) 10, 25,25 

10 R=N 

DO 20 1=1, N 
DO 15 J=1,N 
15 R (1,0) =0.00 
20 R (I , I) =1 . DO 
25 IC=0 
11=0 

T=DABS(A (1,1) ) 

DO 35 1=2, N 

IF (T-DABS (A (1, 1) )) 30, 35, 35 
30 II=I 

T=D A BS (A (I, 1) ) 

35 COHTINOB 

IF (11)110,65,40 
40 IC=IC+1 

IF(tl) 45,55,45 
45 DO 50 J=1 , t! 

S=R (1 J) 

E (1 , J)=R (II, 0) 

50 B (II , 0} =S 
55 DO 60 J=1,H 
S=A (1, J) 

A (1,0) =A (11,0) 

60 A (II ,0) =S 
65 P=A(1,1) 

IF (DABS (P) -EPS) 70 ,70 ,75 
70 IND= 1 
D=0. DO 
GO TO 200 
75 DO 80 0=2, H 
80 A (1 , J) =A (1,0)/P 
J?(K) 85,95,85 
85 DO 90 0=1, H 
90 R(1,0)=R(1,0)/P 
95 DO 170 K=2,N 
KK=K— 1 
T=-1. DO 
DO 105 I=K, H 
DO 98 0=1, K« 

98 A (I,K)=A (I,K)-A (1,0) *A (0,K) 

IF (T-DABS (A (I, K) )) 100, 105, 105 
100 T=DABS (A (I,K)> 

II=I 

105 CONTINUE 

IF (II-K) 110,135,110 
110 IC=IC+1 

IF (H) 115,125,115 
115 DO 120 0=1, » 

S=R ( K , 0 ) 

B (K,0)=B (11,0) 

120 R (II,0)=S 


25 


00000010 
00000020 
00000030 
00000040 
00030050 
00000 060 
00390070 
00000 OA 0 
00000090 
00090 100 

000001 i 0 
00000 120 
00000 130 
00000140 
00000150 
00000160 
00000170 
00000180 
00000190 
00000200 
0000021 0 
00000220 
00000230 
00000240 
00000250 
00000260 
00000270 
00000250 
00000290 
00000300 
0000031 0 
00000 320 
00000330 
00000340 
00000350 
00000360 
00000370 
00000380 
00000390 
00000400 
00000410 
00000420 
00000430 
00000440 
00000450 
00000460 
00000470 
0000048 0 
00000490 
03000500 
00000510 
00000520 
00000530 
00000540 
000005 50 
09000560 
00000570 
00000580 
00000590 
00000600 
00000610 



i 




125 DO 130 J = 1, N 

S=A (K , J) 

A (K, J )=A (XI, J) 

130 A (II,J) = 5 
135 DT=A (K,K) 

IF (DABS (DT) - EPS) 70, 7 0, 140 
140 P=P*DT 

IF(K-N) 1 45,155,145 
145 KP=K+ 1 

DO 150 J=KP,N 
DO 148 1=1, KM 

148 A (K,J) =A (K, J) -A (K,I) *A (I, J) 

150 A (K, J) = A (K, J)/DT 
155 I F ( M) 160,170,160 
160 DO 165 J = 1 , « 

DO 162 1=1, KB 

162 R (K , J)=R (K,J)-A(K,I) *R (I,J) 

165 R (K, J) =R (K, J)/Dr 
170 CONTINUE 

I F (BOD (IC , 2) ) 175,180,175 
175 P=-P 

180 D=P 

IF (H) 185,200,185 
185 II=N 

DO 190 K=2 , N 
KP=I I 

II = II- 1 

DO 190 J = 1 , H 
DO 190 I=KP, K 

190 E (II, J) =B (II,J) - A(II, I) *R (I J) 
200 REIDBU 

END 
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00000 640 
00000 650 
00000660 
00000670 
000006K0 
00000690 
00000700 
00000710 
00000720 
00000730 
00000740 
00000750 
00000760 
00000770 
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00000790 
00000 800 
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Subroutine DCUBIC 


Description : 

Calculates the roots of a cubic equation . 


Argument List: 

Q, R, NRE 

C The coefficients of the cubic equation 

X 3 + C(1)X 2 + C (2 ) X + C (3) = 0 
R Roots, if one real root, it is R{1) 

NRE The number of real roots 


Called by 

. DQRTIC 
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SUBROUT IN E DCUBIC (C,R,NRE) 

00000010 

c 



00000020 

c 



00000 030 

c 


SOLVES POLYNOMIAL EQUATION OF THP TYPE 

00000040 

c 


X**3+C (1) *X* *2 +C (2)*x+ c (3) =0 

00000050 

2 



00000060 

c 


THE COEFFICIENT OP A * *3 IS ASSUMED TO BE 1 

00000 070 

c 



00000080 

c 


R CONTAINS THE: ROOTS 

09000090 

'■' 



00000 100 

c 


NRE CONTAINS THE NUMBER 0? REAL ROOTS 

00000110 

c 



00000 120 

c 


IF THERE IS ONE PEAL ROOT IT WILL BE 

00000 130 

c 


IC P ( 1 ) , WITH THE COMPLEX ROOTS R(2)+-R(3)*I 

00000 140 

c 



00000 150 



DIMENSION C (3) ,R (3) 

00000160 



DOUBLE PRECISION C , R , Cl SQ, P , Q , DEL ,T, A , CRTA ,CRTB , PHI3 , CON, Y 

03000170 



DOUBLE PRECISION B,SQ,HQ 

00000 180 



C 1 SQ=C (1) **2 

00000190 



P = C (2)-C1Sg/3.D0 

00000 200 



Q— C (3) - (C (2) /3.D0-2. D0*C1SQ/27. DO) *C (1) 

00000270 



DEL=4. D0*P**3+27 . 00*0*0 

00000220 



T=C ( 1 )/3 . DO 

00000230 


5 

IF (DEL) 20,10,10 

00000240 


10 

SQ-DSQH? (DEL/1 08. DO) 

00000 250 



HQ=.5D0*Q 

00000200 



A=-HQ+SQ 

00000270 



B=-HQ-SQ 

00000280 



CRTA=DSIGN (CABS (A)**. 33333 33 333333333 3D0, A) 

00000290 



CRIBsDSIGN (DABS (B) 333 333 3 3333 33333 3D 0, B) 

00000300 


15 

Y=CRTA+CRTB 

00000310 



R (1) = Y -T 

00000 320 



R (2) =-.5D0*Y-T 

00000330 



E (3) =. 86 6 025403784 4 3 865 DO* (CRTA -CRTS) 

00000340 



NRE= 1 

00000350 



GO TO 40 

00Q0Q360 


20 

PHI3=BAT AN2 (DSQRT{-D2L/27 .DO) ,-Q)/3.DO 

09000370 


22 

CON=2.00*D3QRT (-P/3- DO) 

00000380 


30 

R (1) = CO N * DCO S (PHI3) -T 

000003 90 



R (2) = -CON*DCOS (1.0471 9755 1 19 6 5977D0-PSI 3) -T 

00000400 



R (3) =-CON*DCOS (1.0471 97551 1965977 DO + PHI 3) -T 

00000410 



NRE= 3 

00000420 


40 

RETUBN 

00000 430 



END 

00000440 
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subroutine DQRTIC 


Description: 

Calculates the roots of a quartic equation. 

Argument List: 

C, R, NRE 

C 

R 

NRE 

Called by: 

. SHADOW 

Calls Subroutines : 

DCUBIC 


The coefficients of the quartic equation 
X 4 + C(1)X 3 + C ( 2 ) X 2 + C ( 3) X + C (4) = 0 
Roots, if 2 real, they are R(l) and R(2) 
The number of real roots 


1 
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SUBROUTINE CQRT IC (C, R, N RE) 


oooooom 

■a 




00000020 

- 




00000030 

c 


SOLVES POLYNOMIAL EQUATION OF TH?? TYPE 


00000040 

c 


X**4 +C (1)*X**3 +C (2)*X**2+C (3) *X+C (4)=0 


00000050 

c 




00000060 

c 


THE COEFFICIENT OF X**4 IS ASSUMED TO BE 1 


000 00070 

3 




00000090 

w 


F CONTLINS THE ROOTS 


00000090 

»*■ 




00000100 

c 


NRE CONTAINS THE NUMBER OF REAL ROOTS 

% 

00000 110 

c 




00000120 

c 


Ii IHFFE ARE TWO REAL ROOTS THEY WILL BE IsT 


OOOOO 130 

c 


R (1) AND R (2) , WITH THE COMPLEX ROOTS R(3)+-R(4)*I 


00000140 

c 




00000150 

c 


IF THERE ARE HO REAL ROOTS, THE COMPLEX 


00000160 

c 


ROOTS ARE R ( 1 ) +—R (2)*IANP R (3J + -R (4) *1 


00000 170 

c 




00000180 



DIMENSION C (4) ,R (4) ,CP(3),Y(3) 


00000190 



DOUBLE PRECISION C, R , CP , Y , C 1 SQ, A ,B, D ,E, F,REAL,D SC R , R AD 


00000 200 



C 1 SQ =C (1) **2 


00000270 



cp 0) =-c (2) 


00000220 



CP(2>=C (1) (3) -4 . D0*C (4) 


00000230 



CP (3) = (4. D0*C (2)-ClSQ) *C (4)-C (3) **2 


00000240 


5 

CALL DCUBIC (CP, Y, DOM) 


00000250 


8 

A=C 1 SQ/4 . DO-C (2)«-Y (1) 


00000260 



B=.5D0*C (1) *Y (1) -C (3) 


00000270 



D=.25D0*Y (1)*Y ( 1 )-C (4) 


00000280 



IF (A) 10, 10, 15 


00000290 


10 

E = 0 , 


00000300 



GO TO 20 


00000310 


15 

E= DS Q ST (A) 


00000320 


20 

IF (D ) 25,25,30 


00000 330 


25 

F=0. 


00000340 



GO TO SO 


00000350 


30 

F=DSIGN (DSQRT{ D>,B) 


00000360 


50 

NRE=0 


00000370 



RE AL=- . 2 5D0+C (1) . 5D0*F 


00000380 



DSCR=RSAL*REAL-.5D0*Y (1)+F 


00000390 


53 

RAD=DSQ3T [DABS ( DSCR) } 


00000400 



I F (DSCR) 60,55,55 


00000410 


55 

NR 3= 2 


00000420 



R (1)=R3AL+RAD 


00 000 430 



R (2)=PEAL-SAD 


00000440 



GO TO 65 


00000450 


60 

R (3) =RSAL 


00000 460 



R (4) = R AD 


00000470 


65 

RE AL=F,EA L~E 


00000480 



DSCR=P.SAL*R.EAL-.5D0*Y (1) -F 


00000490 


63 

RAD=D SQRT (GABS (DSCR) ) 


00000500 



IF (DSCR) 80, 70, 70 


00000510 


70 

N FE= N RE+ 2 


00000520 



R (NRE) =R3AL-RAD 


00000530 



R (NEE- 1)=REAL + RAD 


00000540 



GO TO 90 


00000550 


80 

R ( N RS+ 1 ) =REAL 


00000560 



R (NR2 + 2) =R AD 


00090510 


90 

RETURN 


00000580 



END 


00000590 


30 


3 
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Subroutine EVALMP (EVALMPSS) 

Description : 

Evaluates M and ~ (see Eqs. (3.23) and (3.36)). The form of M 
82" 

which is coded in EVALMP is that shown in Ref. 5, Table 4. It was also 

(i 1*1 

from this form that -g-j- was calculated, while holding F constant when 
taking partials with respect to a, h, k, p, q. 

Argument List: 

X, THETA. MU, AM, EAM, IMFLIG 


X 

Five equinoctial orbital elements 


THETA 

Eccentric longitude 


AMU 

Gravitationsl constant, y 


AM 

3z 

The matrix M = — (Eq. (3.24)) 
3r 


PAM 

The partial of M with respect to the 

orbital elements 

IMFLAG 

Flag, i f = 1, AM and ?AM are calculated; if = 2, only 


M is calculated; if= 3, only PAM is 

calculated 


Common Areas : 

0RBIT2/X1, Yl, RA, PZ20 , PZ26 , PZ29 , PZ35 

Called by: 
rrm 


31 


i 


t 






-1 


i* 

v 

41 


i ; 

M 
i ■ 

j 


J. 



C E VAL M P/I YALM PSS 
C 

C 

C THIS SUBROUTINE EVALUATES THE 5X3 MATRIX M AND THE 
C 5X3X5 PARTIAL OF M WPT X 

r* 

v 

C IF IMFLAG=1, BOTH M (AM) AND ITS PARTIAL (PAH) ARE EVALUATED 

C IF IMFLAG=2 , ONLY H (AM) IS EVALUATED 

C IF IHPLAG=3, ONLY THO, PARTIAL OF M (PAM ) IS EVALUATED 

C . 

c 

SDBBODTTNE SVALMP (X, THETA, AMD, AM, PAM, IHFLAG) 

IMPLICIT REAL*.8 (A-H,0-Z) , INTEGER (I- H) 

DIMENSION X ( 5) , AM(5,3), PAM{5,3,5; 

COMMON /OEBIT2/ Xl , Y1 , RA , PZ20 , PZ26 , PZ29 , PZ35, X 1 DOT, Y 1 DOT 
C 
C 

EN=DSCRT (AMCJ/X (1 )**3) 

RHO= DSQRT ( 1 . DO- T(2)** 2- X( 3)**2) 

BETA= 1 . DO/ (1 . DO +RHO) 

CT= DC OS (THETA) 

ST= DSIN (THETA) 

RA= 1 . Do- X (3)*CP -X (2) *ST 
ZET A= X (3) *ST-X (2)*CT 
BETA 3= B2TA**3/(1 . DO -BETA) 

X1= X(1)*((1.D0 -X (2) **2*3E?A)*CT + X (2) *X (3) *BETA*ST -X(3)) 

Y1= X(1)*((1.D0 -X (3) **2*BETA) *ST ♦ X (2) *X (3) *B2TA*CT -X(2)) 
X11=X1 

yii=ti 

X1DOT= ” { { 1 . DO -1.(2) **2*33VA) *ST -X (2) *X (3) *BET A*CT) *EN*X (1)/RA 
Y 1 DOT- ( ( 1 . DO -X(3)**2*BETA)*CT -X (2) *X (3) *BS TA *ST) *S N* X ( 1) /R A 
PZ1= X (1 ) ^jZETA* (BETA + X ( 2) **2*BETA3) “(X(2)*BETA -ST)*CT/RA) 

PZ2= -X (1 X (-ZETA*X (2)*7, (3) *8ETA3 +1.D0 + (ST -X(2)* SETA) *ST/aA) 
PZ 3= X (1) * (-ZETA*X (2) *X (3) *BETA3-1.D0 +(X(3)*BETA -CT) *CT/RA) 
PZ4= S(1) (-ZETA* (BETA + X (3) **2*BETA3) + (CT -X ( 3 ) * BETA ) * ST/R A ) 
IF (IMFLAG . EQ. 3) GO TO 10 

C IT DO NOT RANT TO EVALUATE PARTIAL OF M , BRANCH TO 10 
AS (1,1) = 2.D0*X1DCT/(EN**2*X (1)) 

AM (1 , 2) = 2.DO*rlDOT/(E»**2*X (1)) 

AH (1 , 3) =0 .DO 
DUM= RHO/(2N*X(l)**2) 

AH (2, 1} = DUM*(PZ2- X (2) *BETA*X 1DOT/ZN) 

AM (2,2) = D UM * ( PZ4 -X ( 2) * B ST A * Y 1 DOT/Z N) 

AM (2,3) = DUM + (X (3) * (X (5) *Y1 -X (4 ) *X 1 ) ) /RHO**2 
A M ( 3 , 1 ) = -BUM* (PZ 1 +X (3) *3ETA *X 1DOT/E N) 

A?! (3,2)= -DUM*(PZ3 + X (3 ) *82T A*Y 1 DOT/EN) 

AM (3 , 3) = -DUM* (X (2) * (X (5) *Y1 -X (4 ) * X 1 ) /RHO**2) 

AM (4, 1) = 0. DO 
AM (4,2) =0. DO 

DUM= ( 1 . DO +X (4) **2 +X(5) **2} / (2. D0*2H*X (1) **2*P.HO) 

AM (4,3) = DUM*Y1 

AH (5, 1) =0. DO 

AH (5, 2) =0. DO 

AH (5, 3). = D U M * X 1 

IP (IMPLAG .EQ. 2) RETURN 

C IF WE ONLY RISH TO EVALUATE! M THEN PROGRAH RETURNS HERE 
10 CA= DSQRT (AMU/X (1))/RA 
PZ5= X (2)*BETA3 
PZ6= X (3)*BITA3 
PZ9= CA*ST/RA 
PZ10= CA+CT/RA 


00000010 
00000020 
00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00000090 
00000100 
00000110 
00000120 
00000130 
00000 140 
00000150 
00000160 
00000170 
00000180 
00000190 
00000200 
00000210 
00000 220 
00000230 
00000240 
00000250 
00000260 
00000270 
00000280 
00000290 
00000300 
00000310 
00000 320 
00000330 
00000340 
00000350 
00000360 
00000 370 
00000380 
00000390 
00000000 
000004 10 
00000420 
00003430 
00000440 
00000450 
00000460 
00000470 
00000480 
00000490 
00000500 
00000510 
00000520 
00000530 
00000540 
00000550 
00000560 
00000570 
00000580 
00000590 
00000 600 
000006 10 
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PZ20= >: (1 ) ii (~2.D0*X (2)*BEIk*CT + X (3)*BE7A*ST + PZ5*ZETA*X(2) ) 
PZ2G= X ( 1 J (X (2} *35TA*.SC -1. DO +PZ6*X (2) *ZETA) 

PZ29= X (1)* (X (3)*BET A*CT -1.90 -PZ5*X (3) *ZE?A) 

PZ35= X (1)* (-2-D0*X ('3j*3ETA*5T + X (2) *BETA*CT -PZ6+X (3) *ZETA) 

PZ11 = -X1D0T/ (2.D0*X (1) ) 

PZ1 2= -Y 1D0T/ (2 . DO*X ( 1 ) ) 

DUM1 = 1 - DO -RA 

P Z 1 3 = -CA* (-2.D0*X (2) *3E£A*ST -X (3)*BE.TA*CT -PZ5*X (2)*DUfl 1 )+PZ9 
1 X1D0T/CA 

PZ 1 4 = -CA* (- X(2)*BE?A*CT -PZ 6 *X (2) *DU M 1 ) +PZ1 0* XI DOT/CA 
PZ15= -CA* (X (3)*.BETA*ST + PZ5*X (3 )*DU K 1) +PZ9*Y1 DOT/CA 
PZ1 6= -CA* (2.D0*X (3)*5ETA*CT +X (2) *BETA*ST + PZ6*DUM1 *X(3) ) 

1 +PZ 10*11 DOT/CA 

DUM= BETA +X (2) *P.Z5 

PZ17= 1.D0+ PZ5*X (2) *( 3. DO/BETA + 1 . DO/ ( 1 .D O-BETA) ) 

PZ18= (2. DO + PZl 7 ) * PZ5 
PZ 1 9= PZ17*PZ6 
D tl M 2 = X{2) *BETA -ST 

PZ21= -X (11* (CT*DUM -ZETA*PZ1 8 4CT*DUM/RA +CT*ST*DUM2/RA**2) 
PZ22= X ( 1 ) (ST*DUM +ZSTA*PZ19 -CT *X (2 ) *PZ6/RA-CT**2 *DUtl2/FA**2) 
PZ23= EETA3* (3. DO/BETA +1.D0/(1.D0 -BETA)) 

PZ24= PZ23*PZ5 
PZ25=PZ23*PZ6 

PZ27= X (1) * {-CT*X (2) *X (3) *SETA3 +ZETA*X (3)* (BETA 3 +X(2)*PZ24) 

1 + (ST| (BETA *X(2) *PZ5 ))/RA +ST**2*DUM2/RA**2) 

PZ2 8= X ( 1 ) (ST*X (2)*X (3)* BETA 3 +ZETA*X (2)" (BETA3 +X (3) *PZ25) 

1 4-x (2)*ST*PZ6/RA +ST*CT*DUM2/RA**2) 

DUM2= X (3 ) *BETA-CT 

PZ30= X (1)* (CT*X (2)*X (3)*BETA3 -ZETA*X(3) (BETA 3 +X(2)*PZ24) 

1 +CT*X^3)*PZ5/RA +CT*ST*DUK2/RA **2) * * 

PZ3 1= X(1) <-ST*X<2) *X (3) BET A3 -ZETA*X(2; (BET A3 +X (3) PZ25) 

1 + CT* ( EETA + X (3)*PZ6) /R A 4-CI**2*DUK2/RA**2) 

DTI 55= BETA +X(3)*PZ6 

PZ32= 1 . DO +PZ6*X (3) * (3. DO/BETA * 1. DO/ (1 . DO. -BETA)) 

PZ 3 3= PZ32*PZ5 

PZ34= PZ32*PZ6 +2.D0*X(3) *3ETA3 

PZ36= X(1)*(CT*DUM -ZETA*PZ33 -ST*X (3)*PZ5/SA -ST**2*DDK2/HA**2) 
PZ37= X (1)* (-ST*DUM -ZET A* PZ 3 4 -ST* (BETA +X (3 ) * PZ6 ) /R A -ST*CT 
1 *DUM2/RA**2) 

DO 20 J= 1 , 2 

20 PAM(1,J,1) = 3.D0*A«(1,J)/(2.D0*X(1)) 

DOM = 2. DO*X (1)**2/AMU 
PAM (1 ,1,2; = PZl 3*DU K 
PAM (1,1,3)= PZl 4* DU M 
RAM (1,2,2)= PZl 5* DU H 
PAM(1,2,3)=PZ16*DUM 
DUM= DSQBT ( AMU*X ( 1 ) ) 

CB=RHO/DOM 

PZ38= — X (2) *CB/RHC**2 
PZ39= -X (3)*CB/RfiO**2 
PAM (2,1,1)= AM (2, 1)/(2.D0*Z(1)) 

BAM (2,1,2)= -CB*BETA*X1 DOT/EH +PZ38* AM (2, 1) /CB + CB* (PZ27 
1 -X (2) *BETA*PZ13/EN -X (2 ) * X 1 DOT* PZ5/EN) 

PAM(2,1,3)= PZ39*AR (2,1) /C3 +.CB*(PZ28 -PZ6*X (2)*'/.! DOT/EN 
1 -X (2j*BSTA*PZ14/EN) 

PAM[ 2,2,1)= AH (2, 2) /(2.D0*X( 1)) 

PAM (2,2,2)= PZ3 8*AM (2 , 2) /CB +CB*(PZ36 -BET A*Y 1D0T/EH -X(2j 
1 *Y1D0T*PZ5/EN -X (2)*BETA*PZ15/EN) 

PAM (2,2, 3) = PZ39*AM (2, 2) /CB +CB*(PZ37 - X (2) *Y1D0T*PZ6/EN 
1 ~X (2)*BETA*PZ16/EN) 

PAM (2,3, 1) = AM (2,3) / (2. DO*X ( 1) ) 
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DfIM1 = X (5) *Y1 -X^4) *X1 

PAM (2 #3,2)- X(3) (X ( 5 )*PZ29 -X (4) *PZ20) / (RHO* DOM) ♦X(2)*X(3) 

1 *DOM1/(RHO**3*DUM) 

PAM (2,3,3) — DnB1/(BHO*DnM) + X (3) *(X (5) *PZ35 -X (4) *PZ26) / (RHO 
1 * DUM ) + X (3)**2*D0N1/(RH0**3*DUM) 

PAM (2,3,4) = -X (3;*Xl/(RHO*DUi1) 

PAH (2,3,5)= X (3) *Y1/(RHO*DriH) 

PAM (3,1, 1 ) = AM (3, 1)/(2.D0*X (1)) 

PAM (3,1,2)= PZ38*AN (3, 1) /CB - CB ^ {PZ2 1 + X (3) *X 1 D OT* PZ5/B H 
1 + X (3)*BETA*PZ13/EN) 

PAH (3,1,3)= PZ39*AH (3, 1) /CB -CB*{PZ22 +(BETA*X1DOT + X(3) 

1 *X1 DOT*PZ6 + X (3 ) *8ETA*PZ14) /EN) 

PAM (3,2,1)= AH (3,2)/(2.D0*X (1) ) 

PAM (3,2,2)= PZ38*AM(3,2)/CB >CB*(PZ30 + X (3)* (Y1 DOT*PZ5 

1 +BETA*PZ15)/ER) 

PAH {3, 2, 3) = PZ39+AM (3, 2)/CB -CB*(PZ31 + (BETA* Y 1 DOT +X (3) 

1 *Y1DOT*PZ6 + X (3) *BETA*PZ16 ) /EN) 

PAM (3,3,1)= AM (3,3)/ (2. D0*X (1) ) 

PAM (3, 3, 2) = -DOM1/ (RHO*DUH) - X (2) * (X (5) +PZ29 -X (4) + PZ20) / 

1 (RBO*DOH) -X (2) **2*DUH1/ (BHO**3*DOM) 

PAM (3,3,3)= -X (2) * (X (5) *PZ35 -X (4) *PZ26) / (BHO*DOM) -X(2)*X(3) 
1 *DGHl/(BHO**3*DDM) 

PAM (3,3,4)= X(2)*X1/(HHO*DUM) 

PAM (3,3,5) = -X (2) *Yl/{RHO+DGM) 

Z5= ( 1 . DO +X(5)**2 *X (4) **2) / (2.DO*DUM*RHO) 

P Z4 0 = -Z5/(2.D0*X (1)) 

PZ41 = X (2) *Z5/RR0**2 
PZ42= X (3) *Z5/BH0**2 
PZ43= X (4)/{DUH*EHO) 

PZ44= X (5)/ (DUM*RHO) 

PAAf(4,3,1)= AM (4, 3 )/ (2. DO *X ( 1) ) 

PAM {4, 3, 2) = PZ41*Y1+Z5*PZ29 
PAM (4,3,3)= PZ42*Y1 +Z5*PZ35 
BAM (4,3,4)= PZ43 + Y1 
PAM (4,3, 5) = PZ 4 4* Y 1 
PAM (5,3,1)- AM(5,3)/(2.D0*X(1)) 

FAM (5,3,2)= PZ4 1 *X1 +Z5*PZ20 
PAM (5, 3, 3) = PZ 42*X1 *Z5*PZ26 
PAM (5,3,4)= PZ43*X1 
PAM (5,3,5) = PZ44*X 1 
DO 30 K= 1 , 5 
PAM(1,3,K)=Q.D0 
DO 30 1=4,5 
DO 30 J=1 ,2 
30 PAH(I,0,iC; =O.DO 
DO 40 1=1,3 
DO 40 J=1 ,2 
DO 40 K= 4 , 5 
40 PAH (1 ,3, K) =0.D0 
RETURN 
END 
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Subroutine FCT (FCTSS) 


, Description: 


Evaluates what is essentially the preaverzged derivative (see 
Section 3.3, Reference 1, Eq. (3.33) , (3.34)) of the state and costate 
due to thrusting. Specifically, 

/SH dt\ = ■ dtt R 2 

VaX dF / = -dF a 


3H dt , „ 3 / dt \ R 

3X dF o\\ dF/J a ( 


a 

The factor — y and the minus sign of Eq . 
R 


(3.34) 


and the factor 


1 

are 

2tr 


taken into account in FUNCT, When FCT is called by OUTP, the thrust 
direction is calculated end then returns prior to the derivative 
calculation. 


Argument List: 


S/ n, 

FI 

F2 

Z 

H 


Z, H, G 

Eccentric longitude 

Another eccentric longitude 

Five orbital elements and their adjcints 

Preaveraged derivative for z. and A corresponding 
to FI 


G 


Preaveraged z and X_ at F ? 


Common Areas: 

A/A, PI 

ORBITl/gECO) , PA 

0RBIT2/X1 , Yl, RA, PR(2,2) , X1D0T , Y1D0T 
. CCOM/ CD44 ) 

SOL/ RS ( 4 ) 

CON/CB, SB, THETA , CTHETA . STHETA . C2T , C4T , S2T , S4T 



Subroutine PCT (FCTSS) (Continued) 


Called by: 

QUAD, FUNCT, OUTP 

Calls Subroutines: 

EVALMP r DRTNI 
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FCT/FCTSS 

THIS SUBPROGRAM IS CALLED BY THE QUADRATURE PROGRAH AND 
EVALUATES THE IKTEGRAHD 


SUBROUTINE FCT (P 1, P2, Z, H, G) 

IMPLICIT REAL*8 (A-H,0-Z) , INTEGER (I-H) 

CCHMON /A/A, AMU, PI 
COMMON /0SBIT1/ YEC(3),PA 

COMMON /0BBIT2/X1 ,X1 , RA , PR (2 , 2) , XI DOT , X 1 DOT 
COMHON /CCOM/CD (4) 

COMMON /SOL/BS (4 ) 

COMMON /COK/CB, SB, THETA ,CT BETA, STHETA ,C2T ,C4T ,S 2T,S4T 
EXTERNAL CORE 

DIMENSION Z (10) ,G (10} ,H{10) , AS (5,3) , PAM (5,3,5) , X (5) ,PRA(5) , 
1 PPA(5) , DRS P (3) , DRS Q (3) , PVEC (3) 


n=o 

F= FI 


EVALUATE M AND PARTIAL OF PI YRT STATE 

DO 5 1=1,5 
X(I)=Z(I) 

0 CALL EV1LHP(X,F,AKU,AH, PAM,1) 

EVALUATE THE COMMON SCALOR FACTOR 

CT=DC0S (F) 

ST=DS IN [P) 

EVALUATE THE PARTIAL OF BA WRT X 
PRA (1) = 0. DO 
PEA (2) = -ST 
PBA (3) = -CT 
PE A (4) = 0 . DO 
PRA (5) * 0.D0 

EVALUATE B TRANSPOSS LAUBDA (PRIMER VECTOR) 

DO 20 1=1,3 
PVEC (I)= O.DO 
DO 20 J- 1,5 

20 PVEC (I)= PVEC CJ>AK (J,I) *Z (J + 5) 

************* 

BRITE (6,1000) .F,Z,CT, ST, EA, PRA ,PYEC 

*********** 

19 A BV EC= O.DO 

DO 22 1=1,3 

22 A BV EC = ABVEC+PVEC (I)* PVEC (I) 

A BV EC= DSQRT ( A BV EC) 

DO 23 1-1,3 

23 PVEC (I> PVEC (iyABV SC 

************* 

WRITE (6, 1000) PYEC , AB VEC 

************** 
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c 

C PRIMER CONE ANGLE. 

.C B = 0 . D 0 
DO 53 1= 1,3 

50 CB= PVEC (I> RS (I>CB 
S8= l.D-10 

IF (DEBS (CB) .LT. 1.D0) SB=DSQRT ( 1 . D 0-CB*CB) 

C 

C PARTIAL OF RS WEST P,Q 

D = 2.D0/(1.D0+Z(4)+Z(4) + Z (5) *2(5)) 

DRSP (1 ) = (-RS(2)*Z(5) -RS (3)) *D 
DRSQ ( 1) = RS (2) *Z (4) *D 

DRSP (2)= RS (1)*Z (5) *D 
D BS Q ( 2 ) = (-RS (1) *Z (4) +RS (3) ) *D 
DRSP (3) = RS(1)*D 

DRSQ (3) = -RS (2) *D 
C 

C THRUST CONE ANGLE 
THG= PI/2. DO 

IF [DABS (SB) .GT. 1.D-8) THG= DATAN2 ( (3.D0*CB+DSQRT (8. DO* SB *5 8 
1 +9.D0*CB*CB) ) ,-4.D0*SB) 

IF (TBG. LT.CD (4) } THG = { 1 . 0000 ID 0*CD (4) ) 

CALL DRTNI {THETA , FT, DPT r CONE ,THG, 1. £-10, 100, IER ) 

TP (IER. GT. 0) URITE (6,1001) IER , T FIG , THET A , FT , D FT 
IP (THETA. GT. PI) WRITE (6,1002) THG, THETA 
IF (THETA. LT.CD (4) ) THETA= CD (4 ) 

STH2T A= DSIN (THETA) 

CTHETA= DCOS (THETA) 

C2T= DCOS (2. DO+THETA) . . 

C4T= DCOS (4 . DO+THETA) 

C 

C 

PCO= STHETA/SB 
RCO = CTHETA-PCO+CB 
DO 7 0 1=1,3 

70 VEC (I)= RCO*BS (I>PCO*PVEC (I) 

A B VEC = 0.D0 
DO 72 1=1,3 

72 A B VEC = ABVEC+VEC (I>VEC (I) 

A BV EC= DSQRT (ABVEC) 

DO 7 4 1=1 ,3 

7 4 VEC (I> VEC (I)/ABVEC 
C 

C ACCELERATION FACTOR AND PARTIAL 

PA= CD ( 1 ) +CD (2) *C2T+CD (3) *C4T 

DUi1MY= 4. DO *C THETA* (CD (2) 4 . DO *CD (3)* C2T) 

DO 80 1= 1,3 
80 PPA (I) = 0. DO 
O.PRSP= 0. DO 
UPRS Q= 0 . DO 
DO 85 1=1 ,3 

DPR S P= UPRSP + VEC (I>DBSP { I) 

85 0 PRS Q= OPBSQ + VEC (I) *DSSQ (I) 

PPA ( 4 ) = D 0 M M Y * U P R S P 
PPA ( 5 ) = DUMMY *UPRSQ 

C 

Q ************* 

C URITE (6,1000) XI ,Y1 , RS ,CB, SB, DRSP, DRSQ, THETA , CTHETA , 

C 1 STHETA,C2T,C4T,S2T,S4T,PA, PPA ,THG, FT, DPT 

C ************ 
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C EV&LTJ&TE FUNCTION 
90 IF (F2.EQ.0.D0) RETURN 

DO TOO 1=1,5 
G(I)= 0. DO 
DO 100 J=1 ,3 

100 G ( I) = G (I) +AH (I, J) *V3C (J) 

HZ = 0. DO 
DO 120 1=1,5 

120 HZ = HZ +Z(I+5)*G(I) 

D0M1 = HZ *RA 

DO 130 1= 1,5 
G (1+ 5) = 0. DO 
DO 130 J=1 ,3 
DO 130 L=1, 5 

1 30**G (1+5) = G (1+5) -Z (L+5) *PAH(L,J,I) *VEC(J) 

Q ** +•********* 

C BRITE (6,1000) G, DUR1,PRA,PPA 

C *’****$******** 

DO 140 1=1,5 

140 G (1+5) = PA* (G (1+5) *RA-HZ*PEA (I)) -PPA (I) *DOM 1 

Q ************** 

C KRITE (6,1000) G,HZ 

C *********** 

C 

DU KMI= R A*P A 

DO 150 1= 1,5 
150 G (I)= G (IJ*D OKMT 
C 

I? (M.EQ.1) RETURN 
DO 200 1= 1,10 
200 H (I) = G (I) 

F=F2 
H= 1 

GO TO IQ 

C1000 FORMAT (1H0, 1P10D12. 4) 

1001 FORMAT (1 HO , 'ERROR IN NERTON ITERATOR FOR COKE ANGLE*/' ISR = 

1 12, 1 THG= , ,1PD20.12,*THETA= • , 1PD20. 12, «FT= »,1PD20.12, 

2 * DFT= * , 1PD20. 1 2) 

1002 FORMAT (1H0, 'ERROR— THETA GREATER THAN 180 DEG'/' THG= *, 

1 1PD20.12, *THETA= ' ,1PD20. 12) 

END 
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Subroutine FUNCT (FUNCTSS) 


i 

! 

t 


Description: 

Called principally by the differential equation integrator to 
evaluate the averaged derivative for the full state and costate. Can 
take into account sun location, shadowing, and oblateness, 

Argument List: 

X, Z, DIR / 

X 

Z 

DERZ 

Common Areas: 

A/A, AMU, PI 

J2/AJ2 

SHAD/FEN, Egg, DFSN (5) , DFEN ( 5 ) , I SHAD 

BCOM/ B (9) 

IS/ ISUN , ISON 

SOL/RSUN ( 3) , Bg 

NU/XNU 

Q/m 


Time 

State and costate 

Derivative of state an d costate 


Called by: 

RK4 , TRAJ 

Calls Subroutines; 

SUN, SHADOW, QUAD, FCT, OBLATE 
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C FUNCT/FUNCTSS 

1 

C SOLAR SAIL 

C THIS SUBROUTINE IS AN INTERFACE BETWEEN THE I !i TEG RAT OR ROUTINE 
C A N 3 THE QUADRATURE: ROUTINE. 

C INCLUDES SHADOW EFFECT. 

2 THIS ROUTINE ADDS THE EFFECT OF OBLATENESS (J2) TO THE DERIV. 

C OBLATE CALCULATES THE EFFECT OF J2. RETURNED AS DZ.J.2, 

C Z IS A VECTOR OF THE AVERAGED STATE AND COSTATZ 
C DERZ IS THE AVERAGED DERIVATIVE OF Z 
C IFCLODE S PENALTY FUNCTION. 

C 

c 

c 

SUBROUTINE F UNCT (X, Z , DERZ) 


IMPLICIT REAL* 8 ( A - H f O-Z ) 


COMMON 

COMMON 




AMU , PI 

J2 


CCHMON /SHAD/ PEN, FEX , DFEN (5),DFEX (5) 
COMMON/BCOH/B (9) 

CCMttON/IS/ISUN, ISON 


COMMON /SOL/RSUH (3) ,RS 
CCHMON /Q/IQ 
COMMON /NU/XNU 


, ISHAD 
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DIMENSION Z (10) , DERZ(IO), G(10), H(10), DZ J 2 ( 1 0) , GEX (1 0) , GEN ( 1 0) 
DIMENSION DZ 1 0 ( 1 0) 

EXTERNAL PCT 

SET UP COEFFS OF COSF AND SINF IN XI AND Y1 AND ? ARTIALS 


BETA- 1.D0/ (1.D0+DSQRT (1-DO-Z (2) **2-Z (3) **2.) ) 
B(l)= 1.D0-Z (2) **2*BZTA 
E (2) = Z (2)*Z (3) *BETA 
B (3) = 1.D0-Z (3)**2*BETA 
EETA3= BETA**3/ (l.DO- BETA) 

A I = Z (2)**2*BETA3 
A2= Z (3)**2*BETA3 
S3= ' B.ETA + A1 
A4= BETA+A2* 

B<4) = -Z (2) (BETA+A3) 

B (5) = Z (3) * A3 
B (6)= -Z (2) *A2 
E (7) = -Z (3) *A1 
B (8)= Z (2) *A4 
B (9) = -Z (3)* (BETA + A4) 


QEX= 0. DO 
QEN= PI+PI 
C 

CALL SUN (X f Z) 
C 
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23 F AC = A / (PI +P I) 

IF {ISON. EQ. 1) FAC=FAC/RS**2 
C 
C 

I S H A D = 0 

I F (ISON. EQ. 0) GO TO 24 
C 

C SHADOW INFLUENCE 

C 

CALL SHADOW (Z) 

. IP (ISHAD.EQ.O) GO TO 24 
QEN = ?EN 
Q EX = FEX 

IF (QEN.LT.QEX) QEN= QEN + 2.D0 + PI 

C 

24 AIF= (QEN-QEX) /DFLOAT (IQ) 

IPV=0 

F1= QEX 
DO 26 1=1,10 

26 DERZ (I)= 0. DO 

27 F2= F1+AIF 

CALL QUAD (FI , F2 , FCT, DZ1 0 , Z, G, H, 1 0 ) 
IF V= IFV + 1 
F 1 = F 2 

DO 37 1=1, 10 

37 DERZ(I)= DERZ{1) + DZ10 (I) 

I F (IFV. LT. IQ) GO TO 27 

C 

31 DO 32 1=1, 10 

32 DERZ (I) = DERZ (I)* FAC 
C PENALTY FUNCTION EFFECT 

IF (XNU.EQ.0.D0) GO TO 38 
E= DSQRTJZ (2)*Z (2) +.2 (3)*Z ( 3)) 

DUH= 1. DO-E 

EL2=XN0/ (Z ( 1 )*Z (l)*DnM*DOM) 

DERZ (6) = -EL2/Z (1) +DERZ (6) 

D0M= EL 2/ (1 .DO—E) 

DOM 1 = 1 . DO 

IF (E.GT. 1.D-10) DOM 1=Z ( 2) /E 
DERZ (7)= DERZ (Vj+DUMI+DDH 
IF (2.GT.1.D-10) D0Ml = Z(3)/2 
DERZ (8) = DERZ (8)+DUM l*DOM 
C 
C 
C 

38 IF (ISHAD.EQ.O) GO TO 90 
C 

C SHADOW INPLUEHCZ 
C 

CALL FCT (QEN , QEX, Z, GEN, GEX) 

JIWX=0. DO 
HWN = 0.D0 
DO 40 1=1,5 
HWX= HHX + Z (1+5) *GEX (I) 

40 HWN= HWN+Z (1 + 5) *GEN (I) 

H WX= HWX*FAC 
HUN= HWN+FAC 
DO 50 1=1,5 

50 DERZ { I + 5)= DERZ (1 + 5) -HHN+DFEN (I) 

C 

C 
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00000850 
00000860 
00000870 
000 00880 
00000990 
00000 900 

00000 91 Q 
00000920 
00000930 
00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 
00001010 
00001020 
00001030 
00001040 
00001050 
00001060 
00001070 

00001 090 
00001090 
00001 100 
00001 1 10 
00001 120 
00001130 
00001140 
00001 150 
00001 160 
00001 170 
00001 180 
00001 190 

HWX*DFEX(I) 00001200 

00001210 

00001220 
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rr **\ 


1 



n o 


90 IF { AJ2. LE. 0. DO) B2TCS8 

OBLATENESS BPPECT 

CALL OBLATE (AJ2 f Z,DZ«J 2,1) 
DO 120 1=1,10 

120 DERZ (I)= DERZ (I)* BZJ2 (I) 
RETURN 
C 

END 


00001 230 
00001 240 
00001 250 
03001260 
00001 270 
00001280 
00001290 
00001 300 
00001 310 
00001 320 



a 



Subroutine INPUT (INPUTSS) 


Description: 

Reads initial data, sets initial values, and prints initial 
information. 

Common Areas: 

XMMM/ ZLO (5) , STEP ( 6 ) , ZERF(6) 

ELEM/ ZPO (5), ZPF (5) 

INT/IPR, IDIM . IDIM2, NIMAX 
TRA/ IFMAX , DT, UEB , EW(10) 

UNITS/UTS, HEM, UIH, UTD , UTKM, DTR , UTKG , UTKW , UTMS2 
T/T&, TO- 
A/A, AMU , PJ. 

WF/ WF (5) 

J2/AJ2 

TC/NOP 

JD/TL. 

IS/ ISUN , ISON 
F/ FLIM , RSTEP 
ORBIT/ NORB 
CONSTR/ICON 
PLOT/ I PL 

TERRA/AE, EC, UU, ENS, AN, ECLMTX (3,3) , EQUMTX( 3,3) 

CCOM/ CD ( 4) 

Q/IQ 

Called by: 

CONTL 

Calls Subroutines : 

PLANET 
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r 

INPOT/INPUTSS 



00000010 






03000020 

W 





00000030 

C 

SOLAR SAIL- - ORBIT RAISING AND ESCAPE 



00000040 

C 

THIS SUBPROGRAM IS CALLED BY CONTI AND READS AND PRINTS 


00000050 

C 

ALL INITIAL DATA AS HELL AS SETS INITIAL CONSTANTS. 


00000 060 

c 

THE UNITS ARE BASED CN INTERNAL KU=1.0, INTERNAL 

DISTANCE 

00000070 

r 

UN 17=1 PLANET RADII. a CIRCULAR 



00000080 

r> 

ORBIT AT 1 PLANET RADII WOULD HAVE A PERIOD OF 2 

PI INTERNAL 

00000090 

-V 

TIME UNITS, 



00000 100 

C 

INPUT 



03000 110 

c 

IPNUH - PLANET NUMBER 1 -M ERCURY , 2- VENDS , 3-EARTH. «- 

MARS 


00000120 

c 


INITIAL ORBIT 



00000130 

c 

A 

(KM) 



00000 140 

c 

V 




00000 150 


I 

(DEG) 



00000 160 

c 

LONG. OF NODE (DEG) 



00000170 

c 

ARGi OF PERICENTER (DEG) 



00000180 

c 


INITIAL GUESSES 



00000190 

c 

LAMBDA A 



00000200 

c 

LAMBDA H 



00000210 

c 

LAMBDA K 



00000220 

c 

LAMBDA P 



00000230 

c 

LAMBDA Q 



00000240 

c 


FINAL CRBIT 



00000250 

c 

A 




00000260 

c 

E 

(NOT USED IF NOP=3) 



00000270 

c 

I 

(NOT USED If NO P= 3 ) 



00000280 

c 

NODE (NOT USED I F NOP=2 OR 3) 



60000290 

c 

.PERICENTER (NOT USED IF NOP=2 OR 3) 



00000300 

c 





00000310 

2 

TF2 (DAYS) , GUBSS FOR FINAL TIRE 



00000 320 

r> 

V 

AC 

(MM/S**2) CHARACTORISTIC ACCELERATION 



00000330 

c 

TL 

JULIAN DATE AT INITIAL TIME 



00000 340 

c 

IRDFLG NC 

IS INAL 


00000350 

c 

1 

END OF INPUT 



00000360 

V 

2 

IPR, PRINT FLAG 

0 


00000370 

c 

3 

NIMAX, MAX. NO. OF ITERATIONS 

20 


00000380 

c 

4 

TFHAX2 (DAYS), MAX. TF 

500. 


00000390 

c 

5 

DT2 (DAYS), TIME STEP FOR D.F. 

10. 


00000400 

c 

6 

0 EB , UPPER ERROR BOUND FOR D.E. 

1-D10 


00000410 

c 

7 

EH, ERROR UEIGHTS FOR D.E. 

1. #1,1 

,1,1,0,... 

00000420 

c 

9 

G M (KK**3/SEC**2) GRAV. CONST. 

SET IN 

PLANETS 

00000440 

c 

10 

NOP = X 5 O.E. SPECIFIED AT TF 

1 


00000450 

c 


= 2, 3 O.E, SPECIFIED AT TF 



00000460 

c 


= 3, SEMIKAJOR AXIS SPECIFIED AT TP 



00000470 

c 

11 

SETS OBLATENESS, A J 2 , = Q.DO 

SET I N 

PLANETS 

00000480 

/*» 

V 

12 

STEP (8) , STEP SIZE FOB NUMERICAL DIFF. 

1,D-6 


00000490 

c 


RSTEP = 0, STEP AS FRACTION IN ITER 

0 


00000500 

c 


= 1, STEP AS CONSTAHT I N ITER 



00000510 

c 

13 

ISON ~ 0, SHACOW EFFECT @FF 

0 


00000520 

c 


= 1, SHADOW EFFECT ON 



00000530 

c 

14 

ISUN = 0 SUM DIST EFFECT CN THRUST OFF 

0 


000005110 

c 


= 1 EFFECT OH 



00000550 

r 

15 

CD { 1 ) , COEFFICIENT IN ACCELERATION FACTOR 

.5 


OOOOO 560 

n 

w 


CD <21 

.5 


00000570 

r 


CD (3 j 

.0 


00000580 

c 

16 

IE, IF = 1, EQUATORIAL FRAME, SOLAR NOTION 

1 


00000550 

c 


= 2, EQUATORIAL FRAME, NO SOLAR MOTION 



00000600 

c 


= 3, ECLIPTIC FRAME, SCLAR ROTION 



00000610 

c 


= 4. ECLIPTIC FRAKE, NO SOLAR MOTION 



00003620 
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/ 




* 5 , 

PLANETAR? FRAME, SOLAR HOTION 




= 6, 

PLANETARY FRANE, NO SOLAR MOTION 


17 

FLIM 

, NORM 

LIMIT IN ITERATION ROUTINE 

1..D-06 

18 

XNU, 

PERICENTER PENALTY FUNCTION FACTOR 

0 . 

19 

NORB 

= 0, NO ORBIT PRINT 

0 



= 1 ,.. 

.,999, ORBIT PRINT ON NORB POINTS 

OF ORBIT 

20 

SETS 

I P L = 1 

SO THAT PLOT DATA IS STORED 

0 

21 

IQ, 

NUMBER 

OF QUADRATURES (1-10) 

2 


SUBROUTINE INPUT 

IMPLICIT REAL*8(A-H,0-Z) , INTEGER (I- N) 

HAMEL 1ST /UN/0TKM,UTS,UTD,UTKG,0TK8,DTHS2 

COMMON /XMMH/ZLO (5), STEP(6) , ZERF (6) 

COMMON /ELEM/ZPO (5) , ZPF (5) 

COMMON /INT/IPR, IDXM, IDIM2, NXSAX 
COMMON /TRA/TFMAX,DT,0E8,E8 (10) 

COMMON /UNITS/UTS ,UTM, OTH, DTD, OTKM ,DTR, UTKG, 0 TK W ,UTMS2 
COMMON /T/TF , TO 
COMMON /A/A, AMO, P I 
COMMON /BF/HF (5) 

COMMON / J 2/ AJ2 
COMMON /TC/NOP 
CCHHON /JD/TL 
COMMON /I S/ISDN, ISON 
COMMON /F/?LIH,K5TEP 
COMMON /ORBIT/ N0S3 
COMMON /CCOM/CD (4) 

COMMON /Q/IQ 
COMMON /PLOT/IPL 

CORHON /TERRA/AE, EC, 00 , ENS , AN , ECLMTX (3,3) ,EQUMTX(3,3) 
COMMON /PLNTS/ IPNIJM 
COMMON /NO/XND 

DIMENSION 8 (5) , PCHAR (4) 

DATA PCHAR/' MERCURY •, 'VENUS ' , 'EABTH ‘,*BARS V 

INTEGER CONSTANTS 

IDIM* 1 0 
IDIM2=5 
I DX M 3 = 6 
I PR=0 
ITF=3 
NIMAX=20 
HOP* 1 
ISON* 0 
ISUN* 0 
KSTEP* 1 
NORB* 0 
IQ* 2 
IS* 1 
I PL* 0 

REAL CONSTANTS 
AMO* 1 . 0D0 


0000063 0 
OOOOO 640 
00000550 
00000660 
00000670 
00000690 
00000690 
00000700 
00000710 
00000720 
ooaoo 730 
00000740 
00000750 
00000760 
00000770 
00000780 
00000790 
00000800 
000008 10 
00000820 
00000830 
00000810 
00000850 

00000 860 
00000870 
00000880 
00000890 
00000900 
00000910 
00000920 
00000930 
00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 
00001010 
00001020 
00001030 
00001040 
00001050 
00001060 
00001070 
00001 080 
00001090 

00001 100 
00001110 
00001 120 
00001 130 
00001 140 
00001 150 
00001 160 
00001170 
00001180 
00001 190 
00001 200 
00001210 
00001220 
00001 230 
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DEB= 1.0D+10 
DO 10 I=1,IDIS2 
EW (I+IDIR2) = 0-D0 
10 3W (I) = I.ODO 

DO 1 2 1= 1 ,IDIR3 
12 STSP(I> 1.D-2 

STEP (IDI R 3} = 1.D-6 
DT2- 10- DO 

DTR= .017453292519943296D0 
PI= 3.141S926535897932D0 
TFMAX2- 500. OEO 
. TFMIN= 0. ODO 
T02=0. ODO 
TL= 0. DO 
FLIM = 1.0-06 
CD (1) = - 5D0 
CD (2) — ,5D0 
CD (3) = 0, DO 
XNU = 0. DO 

ALL READ STATEMENTS FOLLOW 

REAC (5.1006) IPNOM 
URITE (6,2065) PC HA? ( IPNOM) 

IF (IPND11.LT. 1 . OR. IPRtJR. GT.4) GOTO 240 
READ (5.1001) (H (I), I=1,IDIM2) 

WRITS (6,3001) (W (I), 1=1, 1 DIM 2) 

BEAD (5.1001) (ZLO (I), 1= 1 , IDI M2) 

B RITE (6,3001) (ZLO (I) ,1=1,1 DIM 2) 

READ (5,1001) (WF (I),I= 1 ,5) 

WRITS (6,3001) (WF (I), 1=1 ,5) 

READ (5.1001) TF2 

. . -VBIT3 (6,3001) TP2 
READ (5,7001) AC 

WRITE (6,3001) AC 
READ (5.1001) TL 

WRITE (6,3001) TL 

TL MUST BE BETWEEN AEOUT A. E. 1950 AND 2000 
IF ( (TL. LT. 2 . 43 3D 6) .OR. (TL. GT. 2. 452D 6) ) GO TO 230 
CALL PLANETS TO SET UP DEFAULTS 
CALL PLANET (GM) 

>0 BEAD (5,1002) IRDFLG 

WRITE (6,3002) ISDFLG 

IP ( (IRDFLG. GT. 2 1 ) . OR. (IREFLG.1T. 1)) GO TO 200 
>5 GO TO (150,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48, 
1 49,50,51), IRDFLG 

THESE VALUES ARE READ ONLY IP INDICBTEL BY IRDFLG 

32 READ (5,1002) IPB 

WRITS (6,3002) IPR 

IF ( IPR. LT. 0) GO TO 210 
GO TO 20 

33 BEAD (5,1002) NIMAX 

WRITE (6,3002) NIMAX 

IF ( (NIMAX. LT-O) . OR. (BIMAX. GT. 50) ) GOTO 220 
GO TO 20 

34 READ (5,1001) TFHAZ2 

WRITE (6,3001) TFM1X2 

IF ( (TFMAX2. LT.O.DD) .OB. (TPMAX2.GT. 1 ,D3) ) GO TO 220 


00001 240 
00001250 
00001260 
00001 270 
00001280 
00001290 
00001 300 
00001310 
00001320 
00007 330 
00001340 
00001350 
00001 360 
00001 370 
00001 380 
00001390 
00001400 
00001410 
00001420 
00001430 
00001440 
00001450 
00001460 
00001470 
00007 480 
00001490 
00001500 
00001510 
000015’ 10 
00001530 
00001540 
00001550 
00001560 
00001 570 
00001 580 
00001590 
00001600 
00001610 
00001620 
00001630 
00001640 
00001650 
09001660 
00001670 
00001680 
00001690 
00001700 
00001710 
00001720 
00001730 
00001740 
00001 750 
00001760 
00001 770 
00001780 
00001790 
00001800 
00001810 
00001 820 
00007 83 3 
0000i84Q 
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G3 70 20 

00001850 

35 

READ (5,1001) DT2 

00001860 


WRITE (6,3901) DT2 

00001870 

i 

IF < (DT2.LE. 0.E0) .OR. (DT2.GT. 1-03)) GO TO 220 

00001 880 

i 

GC TO 20 

00001890 

1 36 

READ (5,1001) DEB 

00001900 

I 

WRITE (6,3001) DEB 

00001910 

1 

GC TO 20 

00001920 

| 37 

READ (5,1003) EW 

00001930 


WRITE (6,3003) EW 

00001940 


GO TO 20 

00001950 

38 

READ (5.1001) UTKN 

09001960 


WRITS (6,3001) DTK W 

00001970 


GO TO 20 

00001980 

39 

READ (5,1001) GN 

00001990 


URITE (6,3001) GH 

00002000 


GO T 3 20 

00002010 

40 

READ (5,1002) NCP 

03002020 


WRITE (5,3002) NOP 

00002030 


IF { (HOP.LT. 1) .OR. (N0P.GT.3) ) GO TO 220 

00002040 


GC TO 20 

00002050 

4 1 

AG2 = 0. DO 

00002060 


GC TO 20 

00002070 

42 

REAS (5, 1001) (STEP (I) ,I=1,IDIK3) 

00002090 


WRITS (6,3001) (STEP (I) ,I=1,IDIH3) 

00002090 


READ (5,1002) KSTEP 

00002100 


WRITE (6,3002) KSTEP 

00002110 


IF ( (KSTEP. LT. 0) . CR. (KSTEP. GT. 1) ) GO TO 220 

00002120 


GO TO 20 

00002133 

43 

READ (5,1002) ISON 

00002140 


WRITE (6,3002) ISON 

00002 150 


IF ( (ISON. LT.O) .OB. (ISON. GT. 1) ) GO TO 220 

00002160 


GO TO 20 

00002 170 

44 

READ (5,1002) ISUN 

00002180 


WRITE (6,3002) ISUN 

00002 190 


IF ( (ISOH. LT.O) .OR. (ISON. GT. 1)) GO TO 220 

00002200 


GO TO 20 

00002210 

45 

READ (5,1001) (CD Cl), 1= 1 , 3) 

00002220 


WRITE (6,3001 ) (CD(I),I=1 ,3) 

00002230 


IP (DABS (CD (1) +CD (2) +CD (3)-1 . DO) .GT. l.D-8) GO TO 220 

00002240 


GO TO 20 

00002250 

46 

READ (5,1002) IE 

00002260 


YRITE (6,3002) IE 

03002270 


IF ((IS-LT-l).OB. (IE.GT.6) ) GO TO 220 

00002280 


IF (IE.EQ.2 -OR. IE.EQ.4 .OF. IE. EQ.6) ENE=0. DO 

00002290 


I P (IZ.LE.2) GO TO 20 

00002300 


EQOMTX (1 , 1) =1 .DO 

00002310 


EQONTX (2, 1) =0. DO 

00002320 


EQOftTX (3, 1) =0. D 0 

00002330 


EQUMTX (1,2) =0. DO 

00002340 


EQUATE (2 , 2) = 1. DO 

00002350 


EQDMTX (3,2) =0. DO 

00002360 


EQUSTX (I,3)=0.D0 

00002370 


EQDSTX (2 ,3} =0. DO 

00002380 


EQOSTX (3,3) =1. DO 

00002390 


IP (IE.LE.4) GO TO 20 

00002400 


00=0. DO 

00002410 


ECLNTX (1 , 1 ) = 1 . DO 

00002420 


BCLMTX (1,2) =0.D0 

00002430 


ECLKTX (1,3) =0. DO 

00002440 


ECL.Jt.TX (2, 1) =C. DO 

00002450 
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ECLHTX (2 ,2) =1 . DO 
ECLMTX (2 , 3 ) =0. DO 
ECDUTX (3 f 1} =0. DO 
ECLMTX (3,2) =0. DO 
ECLHTX (3,3) = 1, DO 
GO TO 20 

47 READ (5.1003) FLIH 

WRITE (6,3003) FLIK 

IF (FLIM.LT.0.D0) GO TO 220 
GO TO 20 

48 HEAD (5,1001) XNU 

WRITE (6,3001) XNU 
GO TO 20 

49 READ (5,1005) HOBS 

WHITE (6,3005) NORB 

IF { (NOEB.LT.0) ,08. (NOR3. GT. 999) ) GO TO 220 
GO TO 20 

50 X PL= 1 

GO TO 20 

51 READ (5,7002) IQ 

WRITE (6,3002) IQ 

IF ( (IQ.LT.1) .OB. (IQ.GT.10) ) GO TO 220 

GO TO 20 


TIME VALUES ARE CHANGED FBOM CAPS TO OTHER UNITS 

50 UIS= DSQRT(UTKH**3/GH) 

UTK,= DT S/6 0. DO 

UTH=DTS/3600.D0 

UID=0TH/24.DQ 

NOW KOLTIPDY ENE BY DTD, SINCE BE KNOW WHAT IT IS, THIS 
USED TO BE DONE IN EAETH/PLAK2T 
ENE=ENE*OTD 

TO= TO 2/DTD 
TP= TF2/UTD 
TFHAX= TPK1X2/UTD 
DT= DT2/0TD 
T01= TO*OTS 
T F 1 = TF*DTS 
TFHAX 1= TFEAX*UTS 
DT1 = DT*0TS 

KOBE CCflV EBSIONS 

DTHS 2-= OTKH+1.D3/ (CTS**2 
DTKG= 1 » D3 

DTKK= UTKG*0TnS2*UTKM/0TS 

CALCULATE LOWER DIBIT ON THETA 
CD (4 ) =0. DO 

IF ( {CD (2)+CD(3)) .EQ.O.DO) GO TO 155 
OK 3= 8,D0*CD(3) 

CK 1= CD ( 1 ) -CD (2) +CD ( 3) 

CK2= 2- DQ*CD (2)-CK3 . 

IF (CD (3) .EQ.O.DO) CD (4) = DABCOS (-DSQBT (CK 1/CK2) ) 

IF (CD (3) .NE.O.DO) CD(4)= D A R CCS (-DS QRT ( (- 1 . D 0<- DSQRT (1 . DO-4. DO * 


09002 u 60 
00002470 
09002420 
03302490 
00002500 
00002510 
03002520 
00002 530 
00002540 
00002550 
00002560 
00002570 
00002580 
00002590 
00002600 
00002610 
00002620 
00002630 
00002640 
00002650 
00002660 
00002670 
00002660 
00002690 
00002700 
00002710 
00002720 
000027 30 
00002740 
00002750 
00002760 
00002770 
03002780 
00002790 
00002800 
00002 810 
00002820 
00002 830 
00002840 
00002850 
00002860 
00002870 
00002880 
00002890 
00002900 
00002910 
00002920 
00002930 
00002940 
00002950 
00002960 
00002970 
00002980 
00002990 
00003000 
00003010 
00003020 
00003030 
00003040 
00003050 
00003060 
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1 CK1*CK3/CK2**2) ) *CK2/ (2. D0*CK3) ) ) 

155 CONTINUE 


THE PRINTING OF ALL INITIAL VALUES POLLOWS 

WRITE: (6,2000) 

WRITE (6,2065) PC BAR (IPNUfl) 

IF (ISOH.GE. 1) RRITE (6,2042) 

WRITE (6,2002) 

WRITE (6,2003) 

158 WRIT3 (6,2004) (W(I>I=1,5) 

CHANGE FSOa CLASSICAL O.E. TO EQUINOCTIAL O.E. 
DO 160 1=3,5 
160 W (I) = W (I) *DTR 

ZP0{1) = W (1)/ UTKH 

ZPO (2) = W (2)* DSIN(W(5)+W{4)} 

ZPO (3)= W (2)*DC0S (W (5;+W(4)) 

ZPO ( 4 ) = DTAN (W (3) /2.0D0) *DSIN (W (4) ) 

ZPO (5)= DTAN (W (3)/2.0D0) *ECOS (W(4) ) 

A= (AC*1.D-3) /UTBS2 


WRITE (6,2005) 

BRITS (6,2004) (ZPO (I) ,1=1,5) 

BRITS (6,2038) AC, A 
A=A/AE**2 

CD4= (PI- CD (4)) /DTR 

WRITE (6,2063) (CD (I), 1= 1 , 3) , CD4 

IF (ISDN.NE.O) WRITE (6,2044) 

WRITE FINAL CONDITIONS, CHANGE TO EQUINOCTIAL FINAL COND. 

URITE (6,2006) 

ZPF ( 1) = BF (1 )/DTKM 
GO TO (165, 170, 180), NOP 

165 WRITE (6,2003) 

WRITE (6,2004) (WF(l), 1= 1,5) 

DO 165 1=3,5 

166 WF(I) = W F (I) *DTR 

ZPF (2)= BF (2) * DS I N (WP ( , 5)+WF(4)) 

ZPF (3)= WF (2)*DC0S (WP (5) +HP (4)) 

ZPP (4) = DTAN (BF (3) /2.OD0) *DSIN (H? (4)) 

ZPP( 5) = DTAN (BF (3) /2. ODO) *DC0S(WF (4) ) 

DC 167 1=3,5 

167 UP (I>WF (iyDTR 
URITE (6,2005) 

WRITE (6,2004) (ZPF (I) , 1= 1 , 5) 

GO T3 190 

17C ZP? (2) = BF (2) 

ZPF ( 3} = DABS (DTAN (WF (3) *DTR/2. DO) ) 

ZPF (4) = 0. DO 
ZPF (5J=0-D0 
WRITE (5,2031) 

URITE (6,2004) (WF (I) ,1=1 ,3) 


00003070 

00003080 

00003090 

00003100 

00003110 

00003120 

00003130 

000031 40 

00003150 

00003 160 

00003170 

00003180 

00003 190 

00003200 

0000321 0 

00003220 

00003230 

00003 240 

00003250 

00003260 

00003270 

00003280 

00003290 

00003300 

00003310 

00003320 

00003330 

00003340 

00003350 

00003 360 

00003370 

00003380 

00003 390 

00003 400 

00003010 

00003420 

00003 430 

00003440 

00003450 

00003460 

00003470 

00003490 

00003490 

00003500 

00003510 

000 03520 

00003530 

00003540 

00003550 

00003560 

00003570 

00003580 

00003590 

00003600 

00003610 

00003 520 

00003630 

00003640 

00003650 

00 003660 

00003670 
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KBITS 

(6,2032.) 


000036 '33 


WHITE 

(6,2004) 

(ZPF(I) ,1=1,3) 

00003690 


GO TO 190 


00 0 031 0 0 

1 80 

C 3= -Gtt/H F ( 1 ) 


0000371 0 


WRITS 

(6,2064) 

WF ( 1 ) , Z P F ( 1 ) , C 3 

00003720 





00003730 





00003740 

190 

WRITE 

(6,2007) 


00003750 

196 

WGITZ 

(6,201 1) 

ZLO 

00003760 


WRITE 

(6,2008) 


00003770 


FRITE 

(6,2009) 

TF2,TP1 ,TF 

00003780 


WRITE 

(6,2037) 

TL 

00003790 


WRITS 

(6,2036) 

BJ2 

00003800 


WRITS 

(6,2013) 


00003810 


FRITE 

(6,2009) 

o 
fr * 

O 

(N 

O 

00003820 


WRITE 

(6,2015) 


00003830 


KRIT2 

(6,2009/ 

TFMAX2,TF11AX 1 ,TFKAX 

00003840 


WRITE 

(6,2010) 

KSTEP 

000 03850 


WRITE 

(6,2011) 

STEP 

00003860 


PRTTE 

(6,201 E) 


00003870 


WRITE 

(6,2009) 

DT2, DTI , DT 

00003880 


WRITE 

(6,2017) 

DEB 

00003890 


WRITE 

(6,2018) 


00003900 


WRITE 

(6,2019) 

EW 

00003910 


WRITE 

(6, 2058) 

IQ 

00003920 


WRITE 

(6,2020) 

ID IH 

00003930 


WRITE 

(6,2022) 

HIMAX 

00003940 


RRITE 

(6,2048) 

PLIK 

00003 950 


WRITE 

(6,2026) 

ITEM 

00003960 


WRITE 

(6,2027) 

GM 

00003970 


WRITS 

(6,2069) 

XNU 

00003980 


IF (IE 

. EQ.l) WRITE (6,2059) 

00003990 


IF (IE 

.EQ.2) RRITE (6,2060) 

00004000 


IF (IE 

, EQ, 3) WRITE (6,2061) 

000040 10 


IF (IE 

• EQ, 4) WRITE (6,2062) 

00004 020 


IF (IE 

,EQ,5) FRITE (6,2067) 

00004030 


IF (IE 

, EQ, 6) WRITE (6,2068) 

COG04 0 4 0 


WRITE 

(6 ,t)H) 


00004050 


RETURN 



00004060 

200 

WRITE 

(6,2023) 

IRDFLG 

00004070 


STOP 



00004080 

210 

WRITE 

(6,2024) 

IPR 

00004090 


STOP 



00004 1 00 

220 

URITE 

(6,2025) 

IRDPLG 

000041 1 0 


STOP 



00004 120 

230 

#RITE 

(6,2037) 

TL 

00004 130 


STOP 



00004 140 

240 

WRITE 

(6,2066) 

I PROM 

00004150 


STOP 



00004160 

. 




00004170 

1001 

FORMAT 

(F25. 15) 

00004 180 

1002 

FORHAT 

(12) 


00004190 

1003 

FORMAT 

(506.1) 


00004200 

1004 

FORHAT 

(P18.ll) 


00004210 

1005 

FORHAT 

’ (13) 


00004220 

1006 

FORMAT 

(ID 


00004230 

2000 

FORMAT 

‘ (1H1.19X,' MINIK0M TIME OPTIMAL TRAJECTORY 

PROGRAM FOR PLA00004240 

RNETOCENTRIC SOLAR SAIL SPIRAL' / 1H0.50X, • ESCAPE OR CAPTURE') 03004250 

2001 

FORKAT 

' (1H0,53 

X, 1 3H MINIMUM TIME) 

00004 260 

2002 

FORMAT 

* (34H0 THE INITIAL ORBITAL ELEHENTS ARE) 

00004 270 

2003 

FORHAT 

1 (1 HO ,1 0X,6HA (Kfi) ,19X,1KE,20X,7HI (DEG) , 10X, 

16 BLOW A $C NQDE00004280 
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1 (DEG) ,6X, 15H.ARG PERIG (DEG)) 

2004 FORMAT (1P5D23.14) 

2005 PORMAT (1H0,5X,14HA (PLACET RAD) , 16X , 1HH, 22X, 1HK, 22X , 1HP, 22X, 1HQ) 

2006 FORMAT (40H0 THE DESIRED FINAL ORBITAL ELEMENTS ARE) 

2007 FORMAT (32H0 INITIAL GUESSED PARAJETZBS APE) 

2008 FORMAT (21 HO FINAL TIRE ESTIMATE) 

2009 PORMAT (1H , 10X, 1PD22. 1 5, 7H DAYS = , 1 PD22 . 1 5, 1 OH SECONDS -, 1 PD22. 1 


2010 

2011 

2013 

2015 

2016 
2017 
20 18 

2019 

2020 
2022 

2023 

2024 

2025 

2026 
2027 
203 1 
2032 


1 r 6H UNITS) 


FORMAT 

PORMAT 

FORMAT 

FORMAT 

FORMAT 

FOBHAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORHAT 

FORMAT 

FORHAT 

FORHAT 

FORMAT 

FORMAT 

FORMAT 


(50H0 STEP SIZE 
(1P5D23. 14) 


( 1 7H0 
( 1 0 HO 
(27H0 
(36H0 
(35H0 


32H0 
(44HO 
(28H0 
(27R0 
( 1 8 HO 
(11 HO 


FOR NUMERICAL DIPFEREHTIATION, KSTEP = ,I2) 


IS) 


INITIAL TIHE 
TF.MA.X IS) 

TIME STEP ?OR INTEGRATION) 

UPPER ERROR BOUND ON INTEGRATION =, 1PD20.10) 
ERROR' UEIGHTS FOR INTEGB ATION ARE) 

(1P10D8. 1) 

( 1 3 HO DIHEHSION = ,I5) 

MAXIMUM NUMBER OP ITERATIONS =,I5) 

XRDFLG SHOULD BE BETWZEN 1 AND 21, IT IS = ,I5) 
IPR SHOULD BE < 0, IT IS = ,I5) 

BAD INPUT DATA, IRDFLG = ,13) 

1 PLANET RADIUS=,F14.5, 3H KM) 

MO (GH) =,F16.5,13H KH**3/SEC**2) 

(1H0 ,1 0X,6HA (KM) , 19X , 1HE , 20X , 7HI (DEG)) 

(1 HO, 5X, 1 4Rh (PLANET RAD) , 9X , 15HS QRT (H**2+K **2) , 8X , 

1 15HSQRT (P**2 + Q**2) ) 

2037 FORHAT (32H0 JULIAN DATE AT INITIAL TIME IS, PZ0.8) 

2036 FORMAT (6H0 J2 = ,1PD15.7) 

2038 FORMAT (»0 C RAHACTESI3TIC ACCELERATION = • , FI 8. 1 1, 1 HM/S**2 =• 


2042 

2044 

2046 

2048 

2058 

2059 

2060 
2061 
2062 
2063 


1 , 1PD15.7, • 


FORMAT 
FORHAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORYAT 
FORHAT 
FORHAT 
FORMAT 
1 *0 
2064 FORMAT 


(2 4H0 
( *0 
( 28 no 

( 2 3 H 0 
(•0 
{•0 
(•0 
CO 
CO 
CO 


INTERNAL UNITS') 

SHADOW EFFECT IHCLUDED) 

SUN DISTANCE EPFECT INCLUDED') 

INITIAL A (PLANET RADII) =, 1PD25.14) 

NORM LIMIT IN ITER = ,1PD12.5) 

IQ =M3) 

EQUATORIAL FRAME, SOLAR MOTION') 
EQUATORIAL FRAME, NO SOLAR HOTION’) 
ECLIPTIC FRAME, SOLAR HOTION') 

ECLIPTIC FRAME, NC SOLAR MOTION’) 


ACCELERATION FACTOR COEFFICIENTS = 
MAXIMUM THETA (DEG) = « , 1 PD20. 12) 

('0 SEHIHAJOR AXIS= • ,1PD23. 14,* KR = » 


, 1P3D20. 12/ 

1PD23. 14, 


2065 

2066 

2067 

2068 
2069 

3001 

3002 

3003 

3004 

3005 


1 ’ P.R.; 


FORHAT 

FORHAT 

FORMAT 

FORMAT 

FORHAT 

FORMAT 

FORMAT 

FORMAT 

FORXXT 

FORMAT 

END 


C 

c 

CO 
CO 
CO 
{ 1 H 
( 1 H 
( 1 H 
(1H 
( 1 H 


C 3 = 
SPACE 
IPNUH 


1 ,1PD23. 14, ' KM* * 2/SEC** 2 ' ) 
VEHICLE IS IN ORBIT ABOUT » , A 8) 
SHOULD BE BETWEEN 1 AND 4, IT IS 


12 ) 


PLANETARY ORBITAL FRAME, SOLAR MOTION’) 
PLANETARY ORBITAL FRAME, NO SOLAR XOTION’) 
PEEXCENTER PEHALTY FUNCTION PXCTOR =',1PD10.2) 

, F25. 15) 

, 12 ) 

, 5D6. 1) 

, F 1 8 . 1 1 ) 

,13) 


00004290 
00004300 
00004310 
00004320 
00004330 
00004340 
500004350 
00004360 
00004370 
00004380 
00004390 
00004400 
00004410 
00004420 
00004430 
00004440 
00004450 
00004460 
00004470 
00004400 
00004490 
00004500 
00004510 
00004520 
00034530 
00004540 
00004550 
00004560 
00004570 
00004580 
00004590 
00004600 
00004610 
00004620 
00004630 
00004640 
00004650 
00004660 
0000467 0 
00004680 
00004690 
00004700 
000047 10 
00004720 
00004730 
00004740 
00004750 
00004760 
00004770 
00004780 
00004790 
00004800 
00004810 
00004820 
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Subroutine ITER (NRfiS , MODNRSS ) 


Description : 

A Newton- Raphscn iterator. It calculates a nominal trajectory, 
receiving the error in the final condition called It then varies the 

free initial conditions x (including the value of t^) to get a sensitivity 
matrix or partial derivative matrix. By inverting this matrix and pre- 
multiplying new values of x are obtained and a new nominal tra- 
jectory run. This is continued until the norm of the errors in the 
final conditions is less than some inputted tolerance or until the conver- 
gence fails. There is provision for halving the Ax's several times. 

An optional Modified Newton-Raphson iterator does not calculate 
a new partial derivative matrix by running neighboring trajectories at 
each iteration but rather updates the partial derivative matrix using 
the old one and the Lx terms. If convergence is poor a new matrix is 
calculated by running neighboring trajectories, however. 


Argument List: 





KOUNT. NL 

_ FUNCT ■ PRTN 




KOUNT 

Number of trajectory calls 




N I 

The number of iterations 




FUNCT 

Dummy subroutine name " that 
the "y" Cor the iterator (in 

is called 
our case 

to calculate 
this is TRAJ) 

PRTN 

Dummy subroutine name “ used 

t o 

print 

iteration info 


(in our case this is called PRTN) 

Common Areas: 

XKMM/ X ( 5 ) , XS (6) , Y(6) 

INT/IPR, I DIM, IDIM2, MAXNOI 

T/2E., TO 

DY/ PYPT(6) 

F/ KLXM , KSTEP 


Called by: 

CQNTL 

Calls Subroutines: 


DCROUT, PRTN, TRAJ 

i 
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I 


ITER 

NEKTON RAPHSON SBSS 

SPECIAL VERSION FOR 6X6 PARTIAL DZR . MATRIX, 6 DIM. Y. 
SUBROUTINE ITER (KCUNT , NI , FUNCT, FRTN) 

IMPLICIT REA L*8 ( A -H ,0-2) , I NTEGER (I-N) 


X VALUES OF THE INDEPENDENT VAR IA.BLES (I NITIAL, CURRENT, FINAL) 
XS STEP SIZE TO PERTURB X:S TO COMPUTE PARTIAL DERIVATIVES 
Y VALUES OF THE DEPENGENT V ART AELES (CURRENT, FINAL) 
CCMMON/XMMM/X (5) ,XS{6) ,Y(6) 

COMMON /I NT/I PR ,1 TIM , IDI M2, MAXNCI 
COMMON /T/TF , TO 
CCHMON /DY/DYDT (6) 

COMMON /F/FLIH , KSTEP 


C ARRAYS USED INTERNALLY BY THE ITERATOR 

DIMENSION YNOM (6) ,XN (5) ,P (6, 6) ,COEF (6) 

EQUIVALENCE (YNOM, COEF) 

c 

N= 5 
M= 6 
NI = 1 
KCUST=0 
CALL FUNCT 
KCUN!=KOUNT+ 1 
F0=0. DO 
DO 1 5 1=1 ,M 
1 5 FO=FO+Y (I>*2 
10 DO 16 1=1 ,N 
XN (I>X (I) 

16 YNOH (I>Y (I) 

YNOM (M) = Y{M) 

TFN=TF 

CALL PSTN (KOUMT,NI) 

WRITE (6,101 1 )FO 
I F{FO. LS.FLIM) GO TO 90 
IF (NI.GT. MAXNOI) GO TO 80 
C COMPUTE NUMERICAL PARTIAL DERIVATIVES 
WRITE (5,1013) 

DO 17 1=1, K 

17 P (I, M) = DYDT (I) 

DC 25 J=1 , N 
TEMP=X (J) 

STSP=XS (J)*DABS (X (J» 

IF ((DABS (X (J) ) . LT. 1 .D-1 0) .OR . (KSTEP.EQ. 1) ) STEP=XS ( J) 
IF (DABS (X (J) ) .LT. l.D-10) WRITE (6,1014) 

X (J)=X (J)+STEP 
CBLL FUNCT 
WRITE (6,1 000 ) X (J) 

WRITS(6, 1001) (Y (I), I=1,M) 

no 20 1 = 1 ,H 

20 P (I,J)= (Y (I) -YNOM (I)) /STEP 
25 X (J) =TEMP 

KOUNT=KOUNT+N 

WRITE (5,1002) 

DG 30 1 = 1 ,M 


00000010 
00000020 
03003030 
00000040 
00003050 
00000060 
00000070 
00003080 
00000090 
00000100 
00000 1 10 
00000 120 
00000 130 
00000 140 
00000 150 
00000160 
00000 170 
00000180 
00000190 
00000200 
00000210 
00000220 
00000230 
00000240 
00000250 
00000260 
00ODO270 
00001)280 
00000290 
00000300 
00000310 
00000320 
00000330 
00000340 
00000350 
00000360 
00000370 
00000380 
00000 390 
00000400 
00000410 
00000420 
00 000 430 
09000540 
00000450 
00000 460 
00000470 
00000480 
00000490 
00030500 
00000510 
03000520 
00000530 
00000540 
00000 550 
00000560 
00000570 
00000580 
00000590 
00000600 
00000610 
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WRITE (6,1 00 1} (P (I,J) ,J= 1 ,H) 

30 CONTINUE 
DO 35 1=1 ,H 
35 COEF (I)=-YNOtt (I) 

CALL DCROUT (P ,COEF, DET, 0. DO, K , 1 , IND) 

IF (IND.NE.O) GO TO 85 
WRITE (6,1015) DST 
DO 40 1= i , n 

40 IF (DABS (COEF(I) ) .LT. 1.D-12) COEF(I)= 0.D0 
WRITE (6,1003) (COEF (I),X=1 ,N) 

SN= COEF (M) 

WRITE (6,1012) SN 
DO 50 J=1 ,N 
SO X(J) =XN (J) + COEF (J) 

TF=TFN + SN 
IHALV=0 

51 IF (DABS (SN) .GT. (,25D0*TFN) ) GO TO 60 
CALL FUNCT 
KODNT=KOUKT+ 1 
F1=Q. DO 
DO 52 1=1, « 

52 F 1 =F 1 + Y ( I )+ *2 
WRITE (6,1010) FI 
IF(Fl.LT.FO) GO TO 55 

60 YRITE: (6,1006) 

IF (IHALV.EQ. 10) GO TO 95 

IHALV=IHADV+1 

DO 53 J=1,N 

COEP ( J) =COEF (J)/2.D0 

WRITE (6,1000) COEP (J) 

53 X ( J)=XN ( J )+C OEF (3) 

S K= SN/2.D0 
WRITE (6,1012) SR 
IF= TFE + SN 

GO TO 51 

55 IF (NI-HAZKOI) 70,70,80 
70 KI=HI+1 
F0 = F1 
GO TO 10 
80 NI=9999 

WE3TE (6,1006) 

RETURN 
85 NI=9999 

URIT£(6,1007) 

EETDRH 

90 WRITE (6,1005) PO 
RETORH 
95 NI=9999 

RRITE (6, 1009) 

EETURN 

1000 FORMAT (/IX, 1PD23. 15) 

1001 FORM AT (1 X, 1 P5D23. 15) 

1002 PORKAT (21H0PARTIAL DERIV MATRIX) 

1003 FORMAT (1 IHODELXlS ARE/ (1 X ,1 PD23. 1 5) ) 

1005 FORMAT (4H0F0=, 1PD22. 15, 23HCASE CONVERGED. . . FERTIG) 

1006 FORMAT (38H0EXCEEDED MAXIMUM NUMBER OF ITERATIONS) 

1007 FORMAT (16H0MATRIX SINGULAR) 

1008 FORMAT (11H0DELX:S ARE) 

1009 FORMAT (19HOMETHOD CANNOT WORK) 

1010 FORMAT (4H0F1=, 1PD23. 15) 

1011 FORMAT (4H0F0=, 1PD23. 15) 


00000 620 
00000630 
0000064 0 
00000650 
00000 660 
00000 670 
00000680 
00000 690 
09000100 
000007 i n 
00000 720 
00000 730 
00000740 
00000750 
00000760 
00000770 
000007 80 
00000 790 
00000600 
00000810 
00000 820 
00000830 
00000 840 
00000850 
00000860 
00000 870 
00000880 
00000890 
00000 900 
00000910 
00000920 
00000930 
00 000940 
00000950 
00000960 
00000 970 
00000980 
00000990 
00001000 
00001010 
00001020 
00001030 
00001040 
00001050 
00001060 
00001070 
00001080 
00001090 
00001 100 
00001110 
00001 120 
000011 30 
00007 140 
00001150 
00001160 
00001170 
00001180 
00001190 
00001200 
00005210 
00001220 
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1 
I 


1012 FORMAT (IGllO DSL TF = ,1PD23. 15) 

1013 FORMAT (4 OHO X (I>DX (I) FOLLOWED BY CORRESPONDING T) 

1014 FORMAT (2 u HO X (I)= 0, SO DX (I)=XS ( I)) 

1015 FORHAT (15H0 DETERMINER? =,1PD23. 15) 

END 


00001 230 
00001240 
00001250 
00001 260 

00001 270 


i 

(i 
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ITER 

MCDNR/M CDNRSS MODIFIED N - R ITERATOR 


6X6 VSR SION 

SUBROUTINE ITER (KCUNT , NI , FUNCT, PRTN) 

C 

IMPLICIT *EAL*8 (A-H,0-Z) 

C 

C 

» 

C X VALUES OF THE INDEPENDENT VAR TABLES (INITIAL , CURRENT , FINAL) 
C XS STEP SIZE TO PERTURB X:S TO COMPOTE PARTIAL DERIVATIVES 
C Y VALUES OF THE DEPENDENT VARIABLES (CURRENT, FINAL) 
COMMON/XMMK/X (5),XS (6) ,Y (6) 

COMMON /INT/IPR, I DIM , IDIM2 , M AXNOI 
COMMON /T/TF , TO 
COMMON /DY/ DYDT (6) 

CCHHON /F/FL IM , KSTEP 
C 

DIBEPSION YNOK (6), 7. N (5),P (6,6) ,COEF (6),D YDTN (6) 

N=5 
M=6 
18=1 
ICON S=1 
ISN=0 
NI = 1 
KOUNT=0 
CALL PUNCT 
KOONT=KOUNT+1 
F0=0. DO 
DO 15 1=1, M 
15 F0=F0+Y (I)**2 

9 DO 16 1=1, N 

D YDTN (I)= DYDT (I) 

XN (I) =X (I) 

16 YNOH (I) = Y (I) 

YNOM (M) = Y (K) 

TFN=TF 

DYDTK (K) = DYDT (K) 

1 0 CALL PRTN (KOUNT, NI) 
t BRITE (6,1011)F0 

I F(FO.LE.FLIM)GO TO 90 
IF (NI. GT. M AXNOI) GO TO 80 
I? (ISW.NE.O) GO TO 27 
C CCHPUTE NUMERICAL PARTIAL DERIVATIVES 
DO 17 1=1, M 
17 P (I , M) = DYDT (I) 

WRITE (6,1013) 

DO 25 J=1,N 
TEMP=X (J) 

STEP=XS ( j) *DABS (X (J) ) 

IF ( (DABS (X(J) ) ,LT.1.D-10).OR. (KSTEP. EQ. 1) ) STEP=XS (J) 

C IF (DABS (X (J)).XT.1.D-10) WRIT9 (6,1014) 

X (J) = 7 (J)+STEP 
CALL FUNCT 
WRITE (6, 1000) X(*J) 

SHITE (6,1001) (Y (I), 1=1, M) 

DO 20 1=1, M 

20 P (I,J) = (Y (ThYNOM (f)) /STEP 


0000007 0 
00000020 
00000030 
00000040 
00000030 
00000060 
00000070 
00000090 
00000090 
00009100 
00000 1 10 
00000 120 
00000 130 
00000140 
00000150 
00000160 
00000170 
03000 180 
00000190 
00000200 
00000 210 
00000220 
03000 230 
00000240 
00000 250 
00000260 
00 000270 
00000280 
00000290 
00000300 
00000 310 
03000320 
00000330 
00000340 
00000350 
00000 360 
00000370 
OOOOO 380 
00000390 
00000400 
00000 410 
00000420 
00000430 
00 000440 
00000450 
0000046 0 
00000470 
00000480 
000 00490 
00000500 
00000510 
00000520 
00000530 
0000054 0 
00000550 
00000560 
00003570 
00000580 
00000590 
00000600 
000006 10 
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25 

X (J) =TEKP 

00000620 


KCONT=KOONT+N 

00000630 

27 

WRITE (6,1002) 

00000640 


DO 30 1=1, H 

00000 650 


WRITE (6,1001) (P(I,J) , J= 1 , M) 

00000660 

30 

CONTINUE 

000 00670 


DO 35 1=1, H 

00000680 

35 

cost (t;=-ynoh (i> 

00000690 


CALL DCROUT (P, COEP, DET, 0 . DO, M , 1 , IND) 

00000700 


I? (IND-NE.O) GO TO 85 

00000710 


WRITE (6,1015) DET 

00000720 


DO 40 X= 1,11 

00000730 

40 

If (DABS (COEP (I)) .LT. l.D-10) COE?(I) = 0.D0 

00000740 

00000750 


IF (DABS (XN (1) ).5T. 1.D2) GO TO 47 

00000760 


RATS= 1. DO 

00000 770 


DO 45 1= 1,5 

00000780 


RAT= CABS ( COEF(l) ) / {. 8D0+DABS (XN (J^+.IDO) 

00000790 


IF (RAT- GT. RATS) RATS= RAT 

00000800 

45 

CONTINUI 

00000810 


DO 46 1= 1,8 

00000820 

46 

COE? ( I ) = COEP (I) /RATS 

00000830 


WRITE (5,1016) RATS 

00000940 

00000850 

47 

WRITE: (6, 1003) (COEF (I) ,1=1 ,H) 

00000860 


SH= COEP (M) 

00000870 


WRITE (6,1012) SI 

00000880 


DO 50 J= 1 , N 

00000890 

50 

X { J) =SN (J) +COEP {0} 

OOOOOSOO 


TF=TFR +SN 

00000910 


IHALV=0 

00000920' 

51 

CALL FUHCT 

00000930 


K00NT=K0UNT+ 1 

00000940 


F1=0, DO 

00000950 


DO 52 1=1, M 

00000960 

52 

F1=?1+Y (I}**2 

00000970 


WRITS (6,1010) PI 

000 00980 


J?(F1 -LT.FO) GO TO 55 

00000990 


WRITE (6,1008) 

00001 000 


IP(IHALV. EQ. 6) GO TO 95 

00001010 


IHALV=IHALV+1 

00001 020 


DO 53 J= 1 , N 

00001 030 


COEP (J)= COE?(J)/2.D0 

00001 040 


URITE (6,1000) COEP (J) 

00001050 

53 

X (3)=XH (J) +COEF (J) 

00001060 


SH=SN/2. 0D0 

00001 070 


PRITE (6,1012) SH 

00001080 


TP=T?H+SN 

00001090 


GO TO 51 

00001 100 

55 

I F (HI-HAXMOI) 70,70,80 

00001110 

70 

NI=NI+1 

00001120 


IC0NS=RI 

00001 130 


F0=P 1 

00001 140 


SUMDX=0. DO 

00001 150 


DO 76 J=1,H 

00001160 

76 

SUMDX=COEF (J) **2 + SUMDX 

00001 170 


DO 77 1=1, H 

00001 180 


DO 77 J=1,H 

00001 190 


P (I, J) =P (I, J) + (Y (I) *COEF (J) ) /SUHDX 

00001 200 

77 

CONTINUE 

00001 210 


ISW=1 

00001220 
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GO TO 9 
80 N I = 9 9 9 9 

WRITE (6,1006) 

SETtJ EN 
85 KI=9999 

WRITE (6.1007) 

RETURN 

90 WRITS (6, 1005) FO 
RETURN 

95 IF(NI.EQ-1.OR-I8.EQ.10.OR.ICONS-NE.NI)GO TO 100 
IC0NS = IC0N.S t 1 
18*18+1 
DO 96 J= 1,N 
DYDT (J)= DYDIN (J) 

X(J)= XN<J) 

96 Y (J) = YNOM (J) 

Y (K) = YNOK (E) 

DYDT (K) =DY DTN (M) 

TF= TFN 
I S W=0 

URITE (6.1004) 

GO TO 10 
100 111=9999 

BRITE (6, 1009) 

BETURN 

1000 FORHAT (/IX, 1PD23. 15) 

1001 FORM AT(1X, 1P5D23. 15) 

1002 FORMAT (21H0PARTIAL DEBIY KATRIX) 

1003 FORMAT (11H0DELX:S ARE/ (IX ,1 PD23. 15) ) 

1004 FORK AT (35HOFORH HEi! PARTIAL DERIVATIVE MATRIX) 

1005 FORMAT (4H0F0=, 1PD22. 15, 23KCASE CONVERSED. . .FERTIG) 

1 006 FORMAT (38H0EXCEEDED MAXIMUM NUMBER OF ITERATIONS) 

1007 FORMAT (16H0MATRIX SINGULAR) 

1008 FORMAT (11H0DELX:S ARE) 

1009 FORHAT (19H0METF.0D CANNOT YORK) 

1010 FORMAT (4H0F1 = , 1PD23. 15) 

1011 FORMAT (4H0F0=, 1PD23.15) 

1012 FORMAT (10H0 DEL TF =,1PD23.15) 

1013 FORHAT (40HO X(I)+DX(I) FOLLOWED BY CORRESPONDING Y) 

1014 FORMAT (24H0 X(I)=0, SO DX(I)=XS<I)) 

1015 FORKAT (15H0 DETERMINENT =,1PD23.15) 

1016 FORWBT (8H0 RAPS =,1PD23.15) 

END 
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00001400 
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00001480 
00001490 
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00001540 
00001550 
00001560 
00001 570 
00001580 
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00001600 
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00001630 
00001640 
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Subroutine OBLATE 
Description: 

Calculates the oblateness effect on z_ and 
Reference 1) . 

Argument List: 

M2 . Z, DZJ2, U _ 

AJ2 The J 2 coefficient 

Z Orbital elements and adjoints 

DZJ2 The oblateness contribution to 

costate derivative 

IJ Flag, always = 1 

Called by : 

FUNCT 


! 

\ 

1 

i 


{ 

l 

i 


1 

1 

1 

1 
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X (see Section 3.5, 


the state and 





C OBLATE 
C 


C THIS SUBPROGRAM EVALUATES THE AVERAGED DERIVATIVE OF THE 
C STATE AND COSTATE DUE TO THE OBLATENESS, J.2 TERM 
C ASSUMES EARTH EQUATORIAL BADIUS=1, KO=1. IF ROT Cl MUST BE 
C MODIFIED. 

C IF IJ=2, EVALUATE DZ (I) ,1=1,5 ONLY. 

C 


C 

C 

C 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE OBLATE (A J 2 , Z , DZ J2 , IJ) 

IMPLICIT REAL*B(A-H,0-Z) 
DIMENSION Z (10) ,DZJ2 (10) , P J (4,5) 


C1 = 1.5D0*AJ2/Z (1) **3.5 
B1=Z ( 4 ) * * 2 * Z (5) **2 
B7= 1.D0-6. D0*B1 t 3. D0*B1**2 
D2= 1 . DO " Z (2) **2 -Z (3 ) **2 
B2= 1.D0/D2**2 
B6= 1.D0/(1. D0+ E 1) 

C2= B2*B6 
B4= 1 . DO " B1 
c3= C1*C2 
D16= B6*B7 
D3= C3*D1 6 

0252(1)= 0.D0 
DZJ2 (2)= Z (3)*D3 
DZJ2 (3) = -Z (2) *D3 

D4=B4*C3 

DZJ2 (4) = -Z (5) *D4 
DZJ2 ( 5 )= Z (4) *D4 

IF (IJ. EQ. 2) RETURN 

D5= -3. 5D0*C 1/Z (1) 

B8= C2*D16*D5 

PJ (1 ,1)=Z (3) *B8 
PJ(2,1) = -Z (2)*B8 

B 1 2= C2*B4*D5 

PJ (3 , 1) = -Z(5) *B 1 2 
PJ (4,1) = Z (4) *B12 

D2= ,25D0*D2**3 
D7= B6*Z (2)/D2 

B9= C1*D1 6 

PJ (1,2) = Z (3)*B9*D7 

PJ( 2,2)= -B9*(Z(2)*D7 «• C2) 

B 1 3= C1*D7*B4 

PJ (3,2)= -Z(5)*B13 
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00000010 
00000020 
00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00009090 
00000 100 
OOOOOllO 
00000120 
OOOOO 130 
03000140 
00000150 
OOOOO 160 
00000170 
OP 000 1 80 
OG000190 
00000200 
000002 10 
00000220 
00000230 
00000240 
00000250 
00000260 
OOOOO 270 
000002 80 
00000290 
00000300 
00000310 
00000320 
00000330 
00000340 
00000 350 
000003 60 
00000 370 
00000380 
00000390 
00000400 
000004 10 
00000420 
00000430 
00000440 
00000450 
00000460 
00000 470 
00000480 
00000490 
000005 00 
OOOOO 5 10 
00000520 
00000530' 
00000540 
00000550 
00000560 
00000570 
00000580 
030005 90 
00000600 
00000610 
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c 

c 


c 

c 


c 

c 


c 

c 


c 


c 

c 


PJ (4,2)= Z(4)*B13 

D9 = B 6*2 (3)/D2 

PJ (1 ,3) = B9* (2(3) *D9 + C 2) 

PJ (2 ,3) = -Z (2) *B9*D9 

B 14= C 1 *B4*D 9 

P J (3 , 3) = -Z(5)*B14 
PJ (4,3)= Z (4) *B 1 4 

D10= -2.D0*B2*BS**2 

D11=Z (4) *D 1 0 
P 1 2 = -1 2.D0*B4 
P 13= C1*B6 
D 1 5= C2*D 1 2 

B 1 0= D13* (2.D0*B7*D11 + Z(4)*D15) 

PJ (1 , 4) = Z (3) *B 1 0 
PJ (2,4) = -Z (2) *B1 0 

B15= B4*D 1 1 " 2 . D0*Z (4) *C2 

PJ (3,4) = -Z (5) *C 1 *B 1 5 

B 16= C2*B4 

PJ (4,4)= Cl* (B 1 6 r z (4) *315) 

P14= Z (5J*D1 0 
D 1 5= C2*D 1 2 

B1 1= D13* (2. D0*B7*D14 + Z(5)*D15) 

PJ (1,5)= Z ( 3 ) * B 1 1 
PJ (2,5)= -Z (2) *B 1 1 

B 17= P14*B4 “ 2. D0*Z (5) *C2 

PJ(3 , 5) = -Cl* { B 1 6 + Z(5)*B17) 

PJ (4, 5) = Z (4) *C1*E17 


DO 10 J=1 , 5 
DZJ2 (J+5) = 0, DO 
DO 10 1=1,4 

10 DZJ.2 (J + 5) = DZ.J2 (J + 5) -Z (1 + 6) *PJ (I, J) 
BETOBU 
END _ 


00000620 
00000630 
00000640 
00000650 
00000660 
00000670 
00000680 
00000690 
00000700 
00000 710 
00000720 
00 000730 
00000740 
00000750 
00000760 
00000770 
00000780 
00000790 
00000800 
00000810 
00000820 
00000830 
00000 040 
00000850 
00000860 
00000870 
00000880 
00000890 
00000900 
00000910 
00000920 
00000930 
00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 
00002010 
00001020 
00001030 
00001040 
00001050 
00001060 
00001070 
00001080 
00001 090 
00001 100 


62 


1 



Subroutine OUTP (OUTPSS) 


r 


Description: 

Used to print information at each time step of the differential 
equation integration of the trajectory. Optionally it can calculate and 
print information for several points on an individual orbit such as 
thrust direction, sun angles. 


Argument List: 

X, g, DERZ. IHLF, IDIM, PRMT 
T Time 

Z State and costate vector 

DERZ State and costate derivative 

IHLF Flag indicating error in integrator 

IDIM Dimensions of state plus costate 

PRMT Parameter vector containing initial and final time, 

time step, the upper error bound, and a fifth elemsnt 
not used 



Common Areas: 

UNITS /UTS . UTM . UTH , UTD. UTKM . DTR, UT KG. U TKW , UTMS2 
INT/IPR, ID, IDIMZ . NIMAX 
A/^, AMU , P| 

SHAD/FEN, FEX , DFEN ( 5 ) , DFEX(5), I SHAD 

IS/ ISUN , ISON 

SOL/ RS (3) , RSUN 

ORBIT/ NORB 

ORBITl/ UO (3) . PA 

0RBIT2/X1, n, RA, PR(2,2), XIDOT , YIDOT 
CCOM/ CD ( 4 ) 

CON/Cg, SB, THETA . DUB (6) 

PLOT/IPL 


63 


ORIGINAL PAGE Tb 
OF POOR QUALITY 





Subroutine OUTP (OUTPSS) (Continued) 

Called by : 

RK4 

Calls Subroutines: 

SUN, SHADOW , FCT, VCP 




nnnnoono o 


OUTP/ODTPSS 

SCLAB SAIL 

THIS IS TRE OUTP PROG RAM FOR THS R-K 
INTEGRATOR 

equinoctial oh. and costate are used. 

INCLUDES SHADOW TIHE. 


c 


c 


c 


c 


SUBROUTINE OUTP (T ,Z , DERZ , IHLF , ISIS, PR fit) 

IMPLICIT REAL*8(A-H,0-Z)_. INTEGER f I- HI 

COMMON /URITS/OTS ,0TM, UTH, UTD , UTKM , D TR , DTKG , DTK W , DTKS2 
COMMON /INT/IPH,ID,IDIK2 ,NIRAX 

COMMON /A/A, AMU, PI 

COMMON /SHAD/ FEN,FEX # DF3N(5) # DFEX (5) , ISHAD 

COMMON /IS/ISUN, ISOH 

COMMON /SOL/RS (3) ,RSUN 

COMMON /ORBIT/NORB 

COMMON /ORBIT1/U0 (3), PA 

COMMON /ORBIT 2/X1 ,Y1,EA,PE(2,2) ,X1DOT,YlDOT 
COMMON /CCOM /CD (4) 

COMMON /CON /CE, SE, THETA, DUB ( 6) 

COMMON /PLOT/IPL 


DIMENSION PRMT (5) , Z (10) , 
DIMENSION E 1 (3) , E2 ( 3) , 


D5RZ (10) , S (5) 

HI (10) ,G1 (10) , 7 EL (3) 


INTEGER FRSTSS 

DATA PRSTSW /0/ 

IF (IPR.EQ.O) EETURY 

IF (T . EQ. PRKT (1) ) N=0 
51 H+ (y. EQ.PBMT(l)) M=0 


09000010 
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00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00000090 
OOOOO 700 
09000110 
00000120 
00000 130 
00000140 
00000 150 
00000160 
00000 7 70 
00000180 
00000190 
00000200 
OOOOO 210 
00000 220 
OOOOO 230 
00000240 
OOOOO 250 
00000260 
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00000280 
00000290 
00000300 
00000310 
00000320 
00000330 
00000 340 
00000350 


IF ( (T. LT. (. 9999999S99DO*DFLOAT(M) * (FRET (2)-PRYiT (1)) /DFLOAT (IPR) ) 

I „ ,,AND. (IHLF.LT. 11). AND. (T.LT. (. 9999999SSD0*PEBT (2)) ) ) RETURN 
n = M+ 2 

CALL SUN{ T,Z) 

ACR = (A*UTHS2) *1 .D3 
IF (ISOH.EQ.O) GO TO 5 

AC R= ACR/SSUN**2 

TS=0TS*T 

TB=UTM*T 

TH=0TH*T 

TD=UTD*T 

H - 0.D0 

DO 7 1=1 , IDIM2 

H= H + Z (I+IDIH2) *D£RZ (I) 

W (1) = Z(1) *DTK» 

V (2) = 0,DO 

DtJHMY=Z(2)**2+Z (3) **2 

IF (DUMMY. GT. 1. D-40) B (2) =DSQRT (DUMMY) 

H (3) =0. DO 

DDHMY=Z (4 ) **2 + Z(5)**2 

IF (DUMMY. GT. 1. D-40) B(3)= 2.D0*DATAN (DSQRT (DOHKY) ) /DTR 
B (4) = 0.D0 

IF ( (Z (4/.NE.O.DO) .OR. (Z (5) • NE, °* DO) ) « (4) =DATAN2 (Z (4) ,Z(5)) /DTR 

V (5) = 0. DO 


00000360 
00000 370 
00000380 
00000 390 
00000400 
00000410 
00000420 
00000430 
OD 000 440 
00000450 
00000460 
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00000480 
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00000500 
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00000560 
00000570 
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If 





IF ( (7. (2) . HE .O.DO) .OH. (Z (3) . HE. O.DO) ) W (5)=DATAN2 (Z (2) ,Z (3) ) /DTR 
w (5 )= w (5)-w (4) 

IDIil3=IDI3 2+ 1 

C 

C 

C 

C 

WRITE (6,1001) 

YRITE (6,1002) 

WRITE (6,1003) T, TS, TR, TD, N 
WRITE (6,1004) 

URITE (6,1005) (Z (I), 1=1 , IDIH2) 

YRITE (6,1014) 

URITE (6,1005) W 
WRITE (6,1006) 

SHITE (6,1005) (Z (I) , I=IDIM3 , IDIB) 

URITE (6,1007) 

WRITE (6,1005) (DERZ (I), 1=1 , IDIM2) 

URITE (6,1008) 

WRITE (6,1005) (DERZ (I) ,I=IDIH3 , IDIW) 

WRITE (6, 1009) 

PER=2.D0*PI*DSQRT (z (1)**3/AHD) *OTH 
bP= 8 (1) (1.D0+8 (2)) 

PE= 8 (1) * (1. DO-8 (2)) 

C3= (-ASU/Z(1) ) " (CTKM/O-S) **2 
WRITE (6,1010) R, PER, PE, AP, ACR,C3 
C IF PLOT DATA DESIRED, PUT IT OUT 

IF (IPL.NE.O .AND, FRSTS8.EQ.0) WRITE (8) HOBB 
FRST38=1 

I? (IPL.NE.O) YRITE (3) TD,Z , 8, H , PER, PE, AP, ACR, C3 

C 

I P (ISHAD.EQ.O) GO TO 30 

IF (FS2.LT. PER) FEX=£EX*2-D0*PI * 

TSHAD= ( FEZ- PEN +Z (2) (DCOS(PEX) -DCO S (PEN)) -Z (3) (DSIN(PEX) 

1 -DSIN (FEN) ) ) *PER/ (2. DO *PI) 

FPER= TS HAD/PER 
TSH AD= TS HAD*60. DO 
WRITE (6,1013) TSH AD , FPEB 
F EX D= FEZ/DTR 
FEND= FEN/DTR 
PRITE (6,1015) FEND , FEXD 
C 

30 IP (NORB.EQ.O) GO TO 900 
C 

C ********************************** 

c 

WRTT2 (6,1018) RS,RSUN 
50 WRITE (6,1016) 

F1=O.DO 

CP= 2*PI/DFL0AT (NOBB) 

DO 200 1=1 , N0R8 

C FI JS ECC. LONG.: H1,G1 ABE DUMMIES 

CALL ?CT (FI ,0 -DO, Z,H1,G 1) 

C PRIBER CONE AHGLE 

BETA= (PI-DARCOS (CB) ) /DTR 
C THRDST CONE AHGLE 

T RET A D= (PI-THETA) /DTR 
C THRUST CLOCK ANGLE 
VEL ( 1) = X 1 DOT 
VEL (2) = Y1D0T 
VEL( 3)= O.DO 


00000 620 
00000630 
00000640 
00000650 
00000660 
00000670 
00000680 
00000690 
00000700 
00000710 
00000720 
0000073b 
00000740 
00000750 
00000760 
00000770 
00000780 
00000790 
00000800 
00000810 
00000820 
00000830 
00000040 
00000850 
00000860 
00000870 
00000880 
00000890 
00000900 
00000910 
00000920 
00000930 
00000940 
00000950 
00000960 
00000970 
00000980 
00000990 
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00001030 
00001040 
00001050 
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00001070 
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00001090 
00001100 
00001110 
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00001130 
00001140 
00001150 
00001160 
00001 170 
00001 180 
00001190 
00001200 
00001210 
00001220 




U U O U O ft u 


r El XNC E 2 DO HOT HAVE TO BE UNIT VECTORS 
CALL VCP (PS, VEL, E2) 

CALI VCP (E2 ,RS , El ) 

SCAV=0 ,DO 

CCAV=0. DO 

DO 100 J=1,3 

SCAV= SC A V+ U 0 (J) *E2 (J) 

100 CCAV= CCAV+00 (J) *El (J) 

PSIV= DATAN2 (SCAV.CC AV) /DIB 
FD= F1/DTR 

THRUST TO WEIGHT RATIO 

TKR= A*PA*(X1*X 1+Y 1*11) 

WRITE (6,1017) FD,X'1 ,Xl , BETA , THETAD , PSIV, PA,TRE,U0 
F 1- Fl+DF 

IF PLOT DATA DESIRED, PUT IT OUT 

IF (IPL. NE 0) WRITE (8) FD ,X1 ,H £ETA t THETAD, P SIV ,Ph ,TKE,U0 

200 CONTI NOE 

900 RETURN 


00001230 
00001240 
00001 250 
00001260 
00007 270 
00001 280 
00001 290 
00001300 
00001 310 
00001320 
00001 330 
000013*40 
00001350 
00001 360 
00001370 
00001380 
00001 390 
00001 400 
00001 4 10 
00001 420 
00001 430 
00001 440 
00007450 


1001 

1 

1002 


1003 

1004 

1005 

1006 
007 

1 008 
1009 


FORMAT ('0 ************************* ********■**+****♦**************00001 460 
**************************************************************** * J00001470 
FORMAT (5H TIKE, 1 OX , 1 OHTIKE UNITS, 15X,7HSECOHDS, 9X 00001 480 

,5BHOaRS,11X,4EDAYS, 12X,lHN) 00001 490 

FORMAT (1P2D25.7,1P2D15.7,I9/) 00001500 

FORHAT ( « THE EQUINOCTIAL ORBITAL ELERENTS ARE') 00001510 

PORKBT (1P5D22.12/) 00001520 

FORMAT (16HO THE COSTATE IS) 00001530 

FORHAT (32E0 THE DERIVATIVE OF THE STATE IS) 00001540 

FORHAT (34H0 THE DERIVATIVE OF TEE COSTATE IS) 00001 550 

FORMAT (1H0,7X,11 HH A MILTONIAN »8X , 12HPERIOD (HRS),5X, 00001 560 


1 15HPERICENTER (KM) , 7X # 1 4 HAPOCENTER (KK> ,5 X, 1 4HCBAR ACC/RS**2 00601119 

2 , 4X, 1 7HC3 <Kfl**2/SEC**2) ) 

1010 PORK AT (2220**2. IP 4020.1 0 /) 00001590 

1013 FORMAT (24H0 TIRE SPENT IN SHADOW = , F12.6,27H BIN, FRACTION OF P00001 600 
1ESIOD =, F12.8) 00001 610 

( , THE CLASSICAL O.E. •> 00001620 

(1680 EHTRY ANGLE = ,F15.8,14H EXIT ANGLE = ,F15.8) 00001630 

('0 ECC. LONG. XI Y 1 PRIMER 00001640 

THRUST CONE A PSIV ACC FACT T/W RATIO' 00001650 

/50X, *0F»,13X, *UG*,13X, *UW') 00001660 


1014 FORMAT 

1015 FORHAT 

1016 FORMAT 
ICONS A 
2 


1 017 FORHAT (1H0, 1P8D1 5.5/46X, 1P5D 15.5/4 6X, IP 5D 15. 5/ 46X, 1P5D 15.5} 
7018 FORMAT («0 SUN DIRECTION AND CISTANCE: »,1P4D15.5) 

END 


00001 670 
00001630 
00001690 
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Subroutine OUTPC (OUTPCSS) 

Description: 

After the iteration, this subprogram prints a summary of the 
converged iteration parameters, t^, etc. 

Common Areas: 



T/TP, TO 

UNITS/UTS, UTM, UTH, UTD , UTKM. DTR , UTKG, UTKW, UTMS2 
ELEM/ZPO (5) , ZP^yj, 

WF/ WF ( 5) 

TC/NOP 

Called by: 

CCNTL 




f' OtJTPC/OUTPCSS 


_ TBIS SUBPROGRAM WRITES T HE VALUES FOR THE FINAL CONVERGED 
C TRAJECTORY. IT IS CALLED BY TRE MAIN PROGRAH CONTL- 
C SOLAR SAIL- -ORBIT RAISING AND ESCAPE 
C 
C 
c 

SUBROUTINE OUTPC 


C 

IMPLICIT REAL*6 (A-H,0-Z) , INTEGER (I-N) 

C 

COMMON /XMHM/ZLO (5) ,STEP(6),ZERF(6) 

COMMON /Z/ZF (10) ,D2 (10) 

COMMON /T/TF , TO 

COMMON /UNITS/UTS, n ,0TH , UTD , UTKH, DT R, UTKG, UTKK,UTMS2 
COMMON /ELEM/ZPO (5),ZPP (5) 

COMHON / VF/fiF (5) 

COMMON /CCOM/CD (4) 

COIIMON /T C/NOP 

c 

DIMENSION DELZF (5) , DELBF (5) , NFC ( 5 ) 

C 

WFC (1) = ZF ( 1 ) *0TK H 
WFC (2) = 0.0D0 

DD.MMY= ZF (2) **2 + ZF(3)**2 

IF (DUMMY -GT. 1 ,OD-40) WFC (2) =DSQRT (DUMMY) 

WFC (3) = O.QDO 

DUttMY= ZF (4) **2 + ZF (5) **2 

IF (DUMMY. GT. 1 . OD-40) SFC(3)= 2., O.DO*DATAN (DSQRT(DUHHt) ) /DTR 
B.FC (4)- 0.D0 

IF ((DABS (Z F (4) ) . GT. 1 , D-8 ) - AND. (DABS (ZF (5) ) .GT. l.D-6)) 

1 WFC (4) = DATAN2 (ZP (4) ,ZF (5) ) /DTR 
HFC(5)= O.DO 

IF ((DABS (ZF (2) ) . GT. 1 .D-8) .AND. (DABS (ZF (3) ) ,GT. 1 .D-8) ) 

1 RPC (5) =DATAH2 (ZF (2) ,ZF (3))/DTR 

PFC(5)=«?C (5)- 8FC (4) 

DO 10 1=1,5 

DELWF (I) = HFC(I) - SF(I) 

FO DELZ F ( I) = ZF(I) " .ZPF (I) 

TF2= TF*0TD 

TF1= TF*OTS 
C 

RRITE (6,3000) 

RRITB (6,3001) 

URITE (6,3002) WFC 
WFTITE (6,3003) 

WRITE (6,3002) (ZF(I) ,1=1,5) 

RRITE (6,3004) 

WRITE (6,3001) 

IF (NOP. EQ. 1) WRITE (6,3002) DELUF 
IF (NOP.EQ.2) WRITE (6,3002) (DELBF (I) , 1=1, 3) 

IF (NOP.EQ.3) WRITE (6,3002) DELBF (1) 

GO TO (20,30,40) , NOP 
C 

20 YRITE (6,3003) 

WRITE (6,3002) DELZF 

GO TO 100 
C 

30 DELZF (2) = DSQRT (ZF (2) **2+ZF (3) **2) -ZPF (2) 

DELZF (3)= DSQRT (ZP (4) **2+.ZF(5)**2) -ZPF (3) 
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00000430 
00000440 
00000450 
00000460 
00000470 
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"a 



WRITE 

(6,301 1) 


00000620 

.T 

WRITE 

(6,3002) 

(DELS? (I) ,1=1 , 3) 

00000630 


GO TO 

100 


00000640 

■■jU 

C 



000006S0 

4* 

4 0 WRITE 

(6,3003) 


00000660 

si] 

WRITE 

(6,3002) 

DELZF ( 1 ) 

00000670 

*■« : 

C 



00000680 


1 0 0 RRITE 

(6,3006) 


00000690 


YRITE 

(6,3002) 

ZLO 

00000700 

% 

WRITE 

(6,3008) 


00000710 


WRITE 

(6,3009) 

TF2,TF1 ,TF 

00000720 

•i * 

RETURN 



00000730 

I] 

3000 FORMAT 

(35H0 ACTUAL PINBL ORBITAL ELEMENTS ARE) 

00000740 


Hi 

*3-: 


I 

r- 

i'; 

‘4 f 

i 

V; 


3001 


(DEG) , 1 OX, 1 8HL0N ASC 


FORHAT ( 1 H 0 , 1 0 X , 6RA (KH) , 18X, 1H£,20X,7HI 
1 (DZG) ,6X, 1 5HARG PS2IC (DEG)) 

3002 FORHAT < 1 P5D23 . 1 5) 

( 1H0, 5X, 1 4HA (PLBNET RAD) , 1 5X , 1HH , 221, 1H K ,22X , 1 HP , 22X ,1 HQ) 

( 5 1 HO THE ERROR IN THE FINAL O.E. IS (ACTUAL " DESIRED)) 

{ 4 6 HO THE CON VEBG3D INITIAL GUESSED PARAMETERS ARE) 

(29H0 THE MINIMIZED FIHAL TIHE IS) 

(1H ,10X,1PD22. 15, 7H D MS =,D22.15,10H SECONDS =, 1PD22. 15, 600000 820 
1H UNITS) 00000830 

3011 FORMAT (1H0,5X,14HA (PLANET RAD) , 8X , 1 5HSQRT (K**2 +K**2) , 8 X , 00000840 

1 15B5QRT (P**2+Q**2) ) 00000850 

END 00000860 


3003 

3004 

3006 

3008 

3009 


FORMAT 

FORMAT 

FORMAT 

FCRMAT 

FORMAT 


NODEOOOOO 750 
00000760 

00000770 
00000780 
00000790 
00000800 
00000810 
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Subroutine PLANET 


Description : 

Sets various values associated with the planet including orbital 
elements, gravitational constant, radius, oblateness coefficient, rotation 
matrices to transform to ecliptic or equatorial coordinate systems. 

Argument List: 

GM Gravitational constant 

Common Areas: 

JD/2L 

TERRA/ar, EH, W, ENE , AN , ECLMTX (3,3) , EQUMTX(3,3) 

UNITS/UTS, UTM, UTH, UTD, UTRM, TVTg, UTKG, UTKW, UTMS2 

PLNTS/ IPNUM 

J2/ AJ2 


i 



tr 
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C PLAN 'ST 
C 

C THIS SUBPROGRAM SETS ?HE VALUES FOR A PLANET'S 0B3ITAL ELEMENTS 
C AND CALCULATES THE KEAN ANOMALY AT THE INITIAL TIME 
C 0, E. TAKEN ? ROM BATTIN, 1964. EPOCH 1960 JAN. 1.5, JD=2436935. 

C CALCULATES COORDINATE TRANSFORMATION MATRICES AND SETS GRAVITY 
C CONSTANT AND PLANET RADIUS. 

C INPUT 

C I PNUM --PLANET NUMBER 1 -ME SCURF, 2-V3NUS, 3-EART H, 4-EARS 

C TL- - INITIAL TIME, BEGIHNING OF LOW THRUST TBAJ 

2 OUTPUT 

C EARTH'S SEHIHAJOR AXIS 

C IX — 2ARTH'S ECCENTBICITY 

o ^—ARGUMENT OF PZRIH. 

C E HE --ME AN ORBITAL HOTION 

C AN— MEAN ANOMALY AT TL 

C ECLHTX — MATRIX FOB CONVERSION TO EARTH'S ECLIPTIC COORD FRAME 

C EQUMTX--MATRIX FOR CONVERSION TO EQUATORIAL PRAHE 

C IMPLEMENTED FOR EARTH ONLY 

C GM--PLANET GRAVITATIONAL CONSTANT 

C UTKM — RADIUS 

C AJ2--CBLATENESS 

C 
C 

c 

SUBROUTINE PLANET (GM) 

C 

IMPLICIT REAL*8 (A-H,0-Z) 

C 

COMMON /JD/ TL 

COMMON /TERRA/ AS , SC , S , E HE , A N , ECLHTX ( 3 , 3) , EQU MT X (3 ,3 ) 

COMMON /UNITS/ UTS, UTM,UTH,UTD,UTKM,DTR, UTKG, UTKW,UTMS2 
COMMON /PLHTS/ IPNUH 
COMMON /J2/ AJ2 

c 

REAL* 8 SEMI AX (4) ,ECCEN (4) , HLP (4) , HORBH (4) , GRAVC (4) , RAD (4) 
REAL*8 J2 (4) ,MLE (4), ANGOBL (4), LNODE (4;,INCL(4) 

C 

C DATA STATEMENTS 

DATA SSMIAX/0.387099D0,0.7 23 3 22D0, 1.00,1.523691 DO/ 

CATA ECCSN/0.205627D0,6.793D-3, 1. 6726D-2, 9. 3368D-2/ 

DATA MLP/76. 83309D0, 1 3 1 . 0033 IDO , 1 02. 25253D 0, 3 35 . 32269D0/ 

DATA MOB B M/4 . 092339 DO, 1.602131D0, 0. 9 8 56 0 9 DO , 0 . 524O33D0/ 

DATA GB A VC/221 81. 6D0, 324860. IDO, 39 8601 . 2D 0, 42 828 . 4 DO/ 

DATA BAD/2435. DO, 6052. DO, 6378. 16 DO, 3393. 4D0/ 

DATA J2/0. DO, 0. DO, 1 . 0827D- 3, 0. D 0/ 

DATA MLZ/222. 62165 DO, 174. 2943 IDO, 100. 15815D0, 25 8.76729D0/ 

DATA ANGOBL/O.DO, 0. DO, 23. 4500,0. DO/ 

DATA L NO DE/4 7.8571 4D0,76, 31972 DO, 0. DO, Y9.2490 3D0/ 

DATA INCL/7. 00399DO, 3, 39423D0, 0.D0, 1 . 84991D0/ 

C 

C 

C SEKIMAJOB AXIS 

A E = S E M I A X (IPHUH) 

C ECCE.N TRI-CITY 

EC= ECCE N (IPNUM) 

C ARGUMENT Of PERIH. 

W = MLP (IPNUM) -LNODE (IPNUM) 

C MEAN ORBITAL HOTION 
ENE = MOPBtl (IPNUM) 
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GM (KM**3S**2) 

GM=GRAVC (IPNUH) 

: RADIUS (KH) 

UTK M= R AD (IPNUH) 

C 52 

AJ2=J2 (IPNUH) 

C 

C MEAN ANOMALY AT EPOCH 

AN= MLE (IPNUM) -HLP (IPNUM) 

0 MEAN ANOMALY AT TIKE TL 

AN= AN+ENE* (TL- 24 3 69 3 5 . DO) 

AN=A N/360. DO 
AN=?.N-IDINT(AN) 

/•» 

AN= AN*360.D0*DTR 
V = H*D7R 

C ENE SILL BE MULTIPLIED BY UTD IN INPUTSS, SINCE UTD DEPENDS 
C ON U IK M, WHICH CAN BE MODIFIED BY THE INPUT STREAM 
ENE= ENE*DT R 

C 

C SET UP ARRAY FOR USE IN CONVERTING TO EARTH'S ECLIPTIC FRAHE 
C 

DCL=DC03 (LNODE (IPNUM) *DTR) 

DSL=DSIN {LNODE (I PNUM) *DTR) 

DCI=DCOS (INCL(IPNOM) * DC R) 

DS I=DSIN (INCL (IPNUH)* DTE) 

C 

ECLMTX (1,1) -DCL 
ECLMTX (2, 1) =-DCI*DSL 
ECLHTX (3 , 1) = DSI*.DSL 
ECLMTX (1,2) =DSL 
SCLMTX (2,2) =DCI*DCL 
ECLMTX(3,2)=-DSI*DCL 
ECLMTX (1 ^3) = 0 . DO 
ECLMTX (2,3)=DSI 
ECLMTX (3 ,3)=BCI 
C 
C 

C SET UP EQUMTX HR USE I N CONVERSION TO EQUATORIAL FRAME 

C 

C CURRENTLY IMPLEMENTED ONLY FOB EAHTH 
EQUMTX (1 ,1) =1.D0 
EQUMTX (2, 1) =0, DO 
EQUMTX (3, 1 )=Q. DO 
EQUMTX (1,2)=0.D0 
EQUMTX (1,3)=0. DO 
IP (IPNUM. NE. 3) GO TO 5 
DUH= ANGOBL (IPNUH) *DT R 
DCOSD=DCOS (DOM) 

DSIND=D5 IN (DUM) 

EQOMTX (2,2) =DCOSD 
EQUMTX (3 ,2) =-DSIND 
EQUMTX (2, 3) =DSIND 
EQUMTX (3 ,3)=DC0SD 
GO TO 10 
C 

5 EQUMTX (2,2) =1 . DO 
EQOMTX (3,2) =0.D0 
EQUMTX (2,3) =0.D0 
EQUMTX (3, 3) =1. DO 

1 0 CCNTINUB 

C 

RE TURN 
END 


ORIGINAL PAG® ® 1 
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Subroutine PRTN (PRTNSS) 

Description: 

Prints information about the Newton iteration. 

Argument List : 

KOUNT . NOI 

KOUNT The number of tra j ectory'"calls 

NOI The number of iterations 

Common Areas: 

XMMM/ X (5) . XS(6), Y(6) 

T /IE, TO 

Called by : 

ITER 




i 
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PRTK/PBT HSS 

THIS PROGRAM IS CALLED 31 THZ ITERATOR ARE PRISTS 


SUBROUTINE PRTH (MOUNT , SO I) 

IMPLICIT REAL*8(A-H,0-Z) 

CCKKON /XMMfi/X (5) ,XS (6) ,Y (6) 

COMMON /T/TF ,T0 

N = 5 

M=6 

URITE (6,1000) 

KBITS (6,1001) NCI, ROUST 
YRITE (6,1002) 

KBITS (6,1003) (X(J) ,J=1 ,N) 

FRITE (6,1004) 

FRITE (6,1003) (Y(J),J=1,H) 

BRITE (6,1005) TF 

1000 FORMAT (29H0 ITER AD. ‘IRBJECTOBY CALLS) 

1001 FORHAT (1H0,I6,5X,I6) 

1002 FORHAT (2H0X) 

1003 FORMAT (IX, 1P5D23. 15) 

1004 FORMAT (2H0Y) 

1005 FORMAT (5HO TF=, 1PD2 6. 1 6) 

RETURN 

END 
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Subroutine QUAD (QUAD4 , QUAD 8 , QUAD 16 ; QUAD32) 


Description : 

These are 4, 8, 16, and 32 point vector Gaussian quadratures used 

for the averaging. 

Argument List: 

XL, gg, Y, g, G, g, g 

XL 
XU 
FCT 
Y 
2 
G 
H 
N 

Called by : 

FUNCT 

Calls Subroutines; 

FCT 


Lower bound 
Upper bound 

Dummy program name “ called to evaluate integrand 

The resultant integrated function 

Orbital elements and adjoints 

Dummy variable 

Dummy variable 

Dimension 


16 


a 



C QUAD 

C THIS IS A MODIFIED QUADRATURE PROGRAM FOR VECTOR VALUED FUNCTIONS. 
C COMPUTES IHTEGREL OP THE FUNCTION G (OR H> OQSR X FROM XL TO XU. 

C RESULT IS Y. BOUTINE USES A 4 POINT GAUSS QUADRATURE. 

C 

V 

c 

SUBROUTINE QUAD (XL, XU , FCT , Y, Z ,G , H , N) 

C 

IMPLICIT FEAL*8 (A-H,0-Z> s INTEGER (I-N) 

DIMENSION Y (1) ,H (1) , G (1) ,Z (1) 

C 

A= „5D0*{XU«-XL) 

B= XU- XL 

C= ,43056815579702629D0*B 
K= 1 

GO TO 50 
10 DO 20 1=1, N 

20 Y (I)= ,173927R2256872693D0*G (I) 

C— . 1699?O52179242813D0*B 
K=2 

60 TO 50 

30 .DO 40 1=1, N 

40 Y (I)= B* (Y (I) + .32607257743127307D0+G (I)) 

RETURN 

50 CALL FCT (A-C, A+C,Z, H, G) 

DO 60 1=1 , N 
60 G(I)=G(I) + H (I) 

GO TO (10,30) , K 
END 
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QUAD/QUADfi 


00000010 

00000020 

THIS IS A MODIFIED QUADRATURE PROGRAM FOE VECTOR VALUED FUNCTIONS. 

00000030 

COMPUTES TEE INTEGRAL 0? THE FUNCTION G (OR H) 

OVER X FROM XL TO XU. 

00000040 

THE RESULT IS Y. ROUTINE USES A 8 POINT GAUSS 

QUADRATURE. 

00000050 

00000060 

00000070 
OOOOO 080 

SUBROUTINE QUAD (XL, XU , FCT , Y, Z , G , H , N) 


00000090 

00000100 

IMPLICIT REAL*8 (A-H,0-Z) 


00000110 

DIMENSION T (1) ,H (1) ,G (1) ,Z (1) 


00000120 

00300130 

A=.5D0*(XU+XL) 


00000 140 

B-XU-XL 


00000 150 

C=. 4 8 01 4492824076812D0*B 


00000 160 

K=7 


00000170 

GO TO 500 


00000 180 

10 DO 20 1=1, H 


00000 190 

20 Y {I>. 5061426 8 145 188 130 D- 1U (I) 


00000200 

C=. 39833323870681 337D0*B 


00500210 

K=2 


00000220 

GO TO 500 


00000230 

30 DO 40 1=1, N 


00000240 

4 0 r (I ) = Y (I ) +.11 1190517 22 668 724D0*G(I) 


00000250 

C=.26276620495816449D0*B 


00000260 

K=3 


00000270 

GO TO 500 


00000280 

50 DO 60 1=1, N 


00000290 

60 Y fJ>Y (I) +. 15685332293394364D0*G (I) 


00000300 

C=. 9 17 17321 247824 90D-1+B 


000-00 310 

K=4 


00000320 

GO TO 500 


00005330 

70 DO 80 1=1, N 


00000340 

e5 Y fi;=B* (Y (I>. 18134189168918O99D0*G (I)) 


00000350 

RETURN 


00000360 

500 CALL FCT(A-C,A+C,Z,H,G) 


00000370 

DO 510 1=1 ,N 


00000380 

510 G{I)=G(I) + H (I) 


00000390 

GO TO (1 0,30,50,70) ,K 


00000400 

END 


00000410 
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QUAD16 


C THIS IS A fiO DI FI ED QOADBATURE IKTEGHATIOR PBOGBAK FOB 
C VECTOR VALUED FUNCTIONS FO OHE VARIABLE . IT IHTEGRATES 
C G (OR H) OVER X FROX XL TO XU. THE RESULT IS I. 

C EVALUATION IS BT A 16 POINT GBUSS QUADRATURE FORKULA. 

C 

SUBROUTINE QUAD (IL, 10 , PC?, Y, Z ,G , H , H) 

C 

I K PL I CIT R E AL* 8 ( A- H , O-Z ) 

DIMENSION 1(1) ,Z(1) i G ( 1) ,H(1) 

C 

A: . 5D0* (XU+XL) 

B= XU-XL 

C= * 49470046749582497C0+B 
R= 1 

GO TO 500 
10 DO 20 1= 1 , N 

20 Y(I) = .13576229705877047D-1*G{I) 

C= ,4722875115366162900*8 
K= 2 

GO TO 500 
30 DO 40 1=1, H 

40 Y (I) = I (I) + ,31126761969323946D-1*G (I) 

C= -4328 156011 93915 87 D0*B 
K=3 

GO TO 500 
50 DO 60 1=1, R 

60 Y (I ) = I(I)+ .47579255841246392E-1*G (I) 

C= , 377702204 1775O152D0 x B 
K=4 

GO TO 500 
7 0 DO 80 1=1, H 

80 T(I)= Y(I)+ .623 144 85627766936D-1*G (I) 

C= .3089381222013218700*8 
K=5 

GO TO 500 
90 DO 100 1=1, N 

100 X (I)=? (I) +.74797 994 408 283 37D-1*G (I) 

C= . 2290083 88828613 69 D0*B 

K= 6 

GO TO 500 
110 DO 120 1=1, H 

120 Y(I)= Y (I) +■ 8457825969750 127D-1*G (I) 

C= , 14080177538962946D0*B 

K=7 

GO TO 500 
130 DO 140 1=1, H 

140 i (I)= I (I) +.9130170‘752246179D-1*G(I) 

C= .4750625491 8818720D-1*B 
K= 8 

GO TO 500 
150 DO 160 1=1, H 

160 y (I) = B+ (Y (I) +. 947 25305227 5342 5D-1 *G (I) ) 

RETURN 

500 CALL FCT (A-C,A+C,Z,H,G) 

DO 510 1 = 1, S 
510 G (I) = G (I) + B (I) 

GO TO (10, 30, 50, 70, 90, 110, 130, ISO), K 

STOP 

END 
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C QUAD 32 

C 

C THIS IS A MODIFIED QUADRATURE INTEGRATION PROGRAH POR 
C VECTOR VALUED FULCTIONS OF ONE VARIABLE. IT INTEGRATES 
C G (OR H) OVER X FROM XL TO XU. ?HE RESULT IS Y. 

C EVALUATION I S BY A 32 POINT GAUSS QUADRATURE FORMULA, 
r 

SUBROUTINE QUAD ( XL, XU , FCT /Y ,Z,G, H,N) 
r 

IMPLICIT REAL*8 (A-H,0“Z) 

DIMENSION Y(1) , Z (1 ) ,G {1 ) , H (1 ) 

C 

A= - 5D 0* { X0- XL) 

E= XU -XL 

C= -49863193092474078D0*3 
K = 1 

GO TO 500 
10 DO 20 1=1 ,N 

20 Y (I) = . 35093050047350483D-2*G(I) 

C= . 4928O575577263417D0"B 
K=2 

GO TO 500 
30 DO 40 1=1 , N 

40 Y (I) = Y (I) + . 81371 97365452835D“2*G (I) 

C= . 482381 12779375322D0*B 

K = 3 

GO TO 500 
50 DO 60 1=1 ,N 

60 !{!)= y (I) + . 12596032654631 030D-1*G (I) 

C= . 46745303796886984D0*B 

K=4 

GO TO 500 
70 DO 80 1=1 , H 

80 Y (I) = Y (I) + . 171 3693145651 07 17D- 1 *G (I) 

C= . 4481 6057788302606D0*B 
K=5 

GO TO 500 
90 DO 100 1=1, N 

100 Y(IJ=Y{I) + .21417949011113340D-1*G(I) 

C= . 42468380686628499D0*B 
K=6 

GO TO 500 
110 DO 120 1=1 ,N 

120 I (I> Y (I) +. 2549902963 1 1 88088D-1*G (I) 

C= . 39724189798397120D0*B 
K=7 

GO TO 500 
130 DO 140 I- 1 , N 

140 Y(I)=Y(I) + . 29342046739267774D- 1*G (I) 

C= .3660910593701448400*3 

K=8 

GO TO 500 
150 DO 160 1 = 1 , N 

160 Y (I) = Y (I) + .32911111383180923 D“1*G (I) 

C= . 33 1 5221 334651 0760D0*B 
K=3 

GO TO 590 
170 DO 180 1 = 1 ,N 

180 Y (I> Y (I) + . 36172897054424253D-1*G (I) 

C= . 29385787862038116D0+B 
K = 1 0 
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150 

i 00 


21c 

220 


23C 

240 


250 

260 


270 

280 


290 

300 


310 

320 

500 

510 


GC TO 500 
DO 200 1=1, N 

Y (I)~ Y (I) + .39096947B93535153D-1*G (I) 

C- .2534499544661^^1000*8 

K=1 1 

GO TO 503 
DO 220 1 = 1, N 

Y (I> Y (I) + .41 655962713473378D-1*G (I) 

C= -21067563806531767DO B 
K=12 

GO TO 500 
DO 240 1=1, N 

Y(I) = Y(I) + .43926 09650220190 6D-1*G (I) 

C= .165934301 14106332D0*B 
K=1 3 

GG TO 500 
TO 260 1 = 1, » 

Y (I) = Y (I) + . 4559693S347881942D-1*G<I) 

C= . 11964368112606854D0 + E 

K=14 

GO TO 500 
DO 280 1=1, N 

Y (I) = Y (I) + . 46922199540402283D-1*G (I) 

C=. .7223598 079139825D-1*B 

K=1 5 

GO PO 500 
DO 300 1=1 ,N 

Y (I)= Y (If t .4 781936003 96 374 30D-1*G (I) 

C= . 24153832843869158D-1*B 

K= 1 6 

GO TO 500 
DO 320 1=1, S 

Y (I) = B* (Y (I) * .48270044257 36 3900D-1*G(I) ) 
RETURN 
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00000850 
00000860 
00000870 
00000860 
00000890 
00000900 
00000910 
00000920 
00000930 
00000940 
00000950 


CALL FCT(A-C,A+C,Z,H,G) 

DO 510 1=1, H 
G(I)= G(I) + H{I) 

GO TO (10,30,50,70,90,110,130,150,170, 
END 


000 00960 
00000970 
00000980 

190,210,230,250,270,290,310) 00000990 

00001010 
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Subroutine RK4 


Description : 

This is 4 point Runge-Kutta integrator which is just the IBM 
Scientific Subroutine Package version without the accuracy checks. (4) 

Argument List: 

PRMT , Y, DERV , NDIM , IHLF, ggg, OUTP . AUX 

fRMT(5) Initial time measured from launch time, time of 
flight, integration step size, final 2 not used. 

Y State and costate (initial as input, final as output) 

DERY Final derivative 

NDIM Dimension of Y 

IHLF Flag = 12 if TF = TO 

= 13 i f DT (TP " TO) < 0 

FCT Dummy subroutine name “ called to evaluate derivative 

GUT? Dummy output subroutine name 

AUX Dummy array used for intermediate calculation 

Called by : 

TRAJ 


©rich^al page is 

or POOR OUAUTY 
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C RK4 


, THIS IS A RUSGE KUTT A INTEGRATOR OF 4TR ORDER- 
& 

C 

SUBROUTINE RK4 (PR HT , T , DE RY , N DIM , IHLF, FCT , 0 U TP , A 0 X) 

C 

IMPLICIT BEAL*8(A-H,0-Z) 

DINENSION Y (1) , DER Y (1) ,AUX (1 f 1 ) , A (4 ) , B (4 ) , C (4) , PRMT ( 1) 
C 

X=PRMT (1 ) 

XEND=PRMT (2) 

H=PRMT (3) 

PRMT (5) =0 . DO 

IREC=0 

IHLF=0 

CALL FCT (X, Y , DERY) 

C 

C ERROR TEST 

IF (H* (XEHD-X) ) 38,37,2 
C 

C PREPABATIOB FOR R-K METHOD 

2 A(l)=,SDO 

A (.2) = . 29289 32 188 134 52 48D0 
A ('3j = 1.7071067811865475DO 
A (4) = , 1666666666 6 666667D0 • 

B (1) =2. DO 
B (2)= 1.D0 
B (3) = 1 -D 0 
E (4) =2. DO 
C (1) =. 5D 0 
C (2} = A (2) 

C (3) =A (3) 

C (4) = , 5D0 
C 

C PREPARATION OF FIRST R-K STEP 
DO 3 1=1 , HDIM 

3 AUX (1 , I)=0. DO 

IEHD=0 

C 

C 

C START OF A R-K STEP 

4 IS=0 

IF ( (X+2.D0+H-X2ND) *H) 7,6,5 

5 H=(XEND-X)/2.D0 

6 IEND= 1 

7 CALL OUTP(X ,Y,DERY ,IREC ^KDIM, PRKT) . 

IF (PRMT (5)) 40, 8, 40 

C 

C 

C START OF R-K LOOP 

8 0= 1 

10 AJ=A (J) 

BJ=B (J) 

CJ=C (J) 

DO 11 1=1 , HDIM 
R1=H*DERY (I) 

R2=AJ*(R1~BJ*A0X (1 ,1)) 

Y (I>Y (I)+R2 
R2=R2+R2+R2 

11 AOX (1,I)=AUX (!,!)♦ R2-CJ*R1 
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OOOOOO 1C 
00000020 
00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00000090 
00000 100 
ooooo no 
00000120 
00000130 
00000 1 40 
00000 150 
00000160 
00000170 
00000 180 
00000190 
00000200 
00000210 
00000220 
00000230 
00000240 
00000250 
00000260 
00000 270 
00000230 
00000290 
00000300 
00000310 
00000320 
00000330 
00000340 
00000350 
00000360 
00000370 
00000380 
00000390 
00000 400 
00000410 
00000420 
00000430 
00000440 
00000450 
00000460 
00000470 
00 000480 
00000490 
00000500 
00000510 
00000520 
00000530 
00000540 
00000550 
00000560 
00000570 
00000 580 
00000590 
00000600 
00000 610 
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Subroutine SHADOW (SHADOWSS) 


Description: 

Sets up the shadow quartic equation (Appendix D, Reference 1)/ 

solves it, checks the roots to find the entrance and exit of the orbit 

clF for 

from shadow if there is intersection and calculates T (Eq. D. 5) 

dz. 

entry and exit. 

Argument List: 

Z. 

Z Orbital elements and their adjoints 


Common Areas: 

SOL/ XSm , vsttm 7.STTN. RSUN 

SHAD/EEM-, E , DFEN ( 5 ) , DFEX ( 5 ) , ISHAD 


Called by: 

FUNCT, OUTP 


PAGE BLANK NOT FTLWEO 
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nnonnnnnn o n n n nnnnoooooononoouoono 


SHADOH/SHADOWSS 

THIS PROGRAY D'ST ERMINES IF h GIVEN ORBIT PASSES THROUGH 
THE EARTH'S SHADOW AND IF SO WHAT THE ENTRY AND EXIT 
ANGLES ARP. IT ALSO CALCULATES THE PARTIAL DERIVITIVE 
VECTOR OP F WRT THE O.E. AT THE EHTRY AND EXIT POINTS, 

INPUT AND OUTPUT IK EQUINOCTIAL COOED, 

INPUT 

Z-~ 10 VECTOR OF O.E. AND COSTATE (COSTATE NOT USED) 

XSUN , YSON ,ZSUN--S0N' S DIRECTION IN EQ. COORD. (UNIT VECTOR) 
OUTPUT 

FEN” -ENTRY AHGLE 
PEX--EXIT ANGLE 

DFEN--D2RI VATI V2 OF P AT EHTRY 

DFEX--DER IVATIVE OF F AT EXIT 

IS HAD--FL AG=0 IF ORBIT NOT INTERSECT SHADOW 

=2 IF ORBIT ENTER AND EXIT FROM SHADOV 


SDBBOUTINE SHADOW (Z) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON /SOL/ ISON , ISDN f ZSON,RSUN 

COMMON /SHAD/ FEN ,FEX,DFEN (5) , DFEX (5) , IS HAD 

DIMENSION DSDX (5),AP(4) ,B?{4) ,Z (10) 

NAMELIST /DUHPEN/ FEN , DOM 
NAMELIST /DUMPEX/ FEX ,DUH 
NAMELIST /EQ/ AP, BT, NRE 

NAMELIST /EQ2/BET A ,B1 , B2 , B3 , D1 , D2, D3, HI , H2, H3, G 1 , G2, AO 
NAMELIST /PR4/ DU M , DS DF , DSDX 
NAMELIST /PR 5/ I , II , CP, SF ,X 1 , 1 1 
NAMELIST /PR6/ I,II,EQN 


CALCULATE POLYNOMIAL COEF. 

BETA = DSQBT (1.D0-Z (2)**2-Z (3)**2) 

BETA= 1 . DO/ (1-DO+BETA) 

B 1 = 1.D0-Z (2)**2*BBTA 

B2= Z (2)*Z (3)*BETA 

B 3= 1.DO-Z (3) **2*BETA 

D 1 = 1. D0-XS0N**2 

D 2= 1 . D0-ESUH**2 

D3 = 2.D0*XSUN*YSUN 

C1=B2**2 

C2 = B 3**2 

C3=B2*B3 

C4=B 1*B2 

H1=D1* <B1 **2-Cl) +D2* (C1-C2) -D3* (C4-C3) 
H2=-2.D0*{D1*B1*Z (3) +D2*B2*Z (2)) +D3*(B2*Z (3)+Bl*Z (2)) 
H3=D1 * (Cl + Z (3) **2) *D2* (C2*Z (2) **2) -D3*(C3+Z (2) *Z (3)) 

1 -1.D0/Z(1 y*2 

G 1= 2. DO* (D1*C4 + D2*C3) -D3* (C1*B1*B3) 

G2= -2. DO* (D1*B2*Z(3)*D2*B3*Z (2))* D3* (B3*Z (3) +B.2*Z (2)) 
Cl= G 1 ** 2 
C2= G 2 **2 

C3= G 1*G2 


00000 010 
00000020 
00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00000090 
00000100 
ooooo no 
00000 120 
00000130 
00000140 
00000 150 
00000 160 
00000 170 
00000180 
OOOOO 190 
00000200 
00000210 
00000220 
00000230 
00000240 
00000250 
00000260 
00000270 
00000280 
00000290 
00000300 
00000310 
00000320 
00000330 
00000340 
00000 350 
00000360 
00000370 
00000380 
00000390 
00000400 
00000410 
00000 420 
00000430 
00000440 
00000450 
00000460 
00000470 
00000400 
00000440 
00000500 
00000510 
00000520 
00000530 
00000540 
00000550 
00000560 
00000570 
00000580 
00000590 
00000600 
00000610 
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AO: Hl**2+Cl 
C WRITE (6,E Q2) 

AP (1)=2. DO* (K1*112+C3) /AO 
AP<2) = (H2**2+2.D0*H3*H1-Cl+C2)/A0 
A P { 3 ) - 2. DO* (H3*H2-C3)/AO 
AP{ 4) = (H3**2-C2) /AO 
C 

C CALL SUBROUTINE TO SOLVE A QUARTIC EQN. 

C 

CALL DQBTIC(AP,RT,NRE) 

C 

C WRITE <6,EQ) 

C 

C NBE= NUMBER OF REAL ROOTS, MUST BE EQUAL TO 0,2, OB 4 
C BOOTS ARE RT (I> OR RT ( 1) ,ET(2) ,RT(3)+-RT (4)*T r OR 
C RT(1) + -RT ( , 2;*I,RT{3}+-fiT<4) *1 

IF ( (NRE. EQ. 1), OR, (NRE.FQ.3)) GOTO 130 

C 

FEN= 100. DO 
?EX= 100. DO 
ISKADs 0 
1=0 

10 1=1+1 

IF ( (ISHAD.EQ.2) ,OB. (I.EQ. (NRE + 1)}} GO TO 120 
C 

CP= RT (I) 

$F= DSQBT (1.D0-CP**2) 

C 

C HEMISPHERE CHECK 

11=1 

20 X1= B1*CF«-B2*SF-Z (3) 

Y 1= B3*5F*B2*CF~Z (2) 

C WRITE (6, PR5) 

IF {{X1*X50N*Y1*YSUN) .LT.O.DO) GO TO 40 
30 IF (II. EQ. 2} GO TO 10 
11=2 
SF=-S? 

GO TO 20 
C 

C IS SHADOW EQUATION ZERO? 

40 EQN= bl*X1**2+D2*Yl**2-D3*Xl*Y1-Z (1) ** (-2) 

C RRLTE <6,PR6) 

I? (DABS (EQH) .GT. 1.D-6) GOTO 30 
C 

C ROOT HAS PASSED TESTS— NOW CHECK TO SEE IF EXIT OE ENTRY ANGLE 
DXDF= -B1*SF+B2*CF 
D YDF= -B2*SF+B3*CF 

DSDF= (2. D0*D1*X1-D3*Y1) *DXDF + (2.D 0*D2* Y1-D3*X1 ) *D YDF 
D0M= DATATI2 (SP, CF) 

IF (DSDP) 70,50,60 
C 

C OBBIT IS TANGENT TO SHADOW 

50 BRITE (6,1010) 

GO TO 30 
C 

C IS PEX ALREADY FOUND? 

60 IF (FEX. EQ. 1.D2) GO TO 80 

C YES 

WRITE {6, CUM FEX) 

GO TO 30 

c 


00000620 
00000630 
00000540 
00000650 
00000660 
00 000 67 0 
00000680 
000 00690 
00000700 
00000710 
00000720 
00000730 
00000740 
00000750 
00000760 
00000770 
00000 780 
00000790 
00000800 
00000810 
00000820 
00000830 
00000840 
00000 850 
00000860 
00000870 
00000880 
00000890 
00000900 
00000910 
00000920 
00000 930 
00000940 
00000950 
00000960 

000 00970 
000009S0 
00000990 
00001000 
00001010 
00001020 
00001030 
00001040 

00001 050 
00001 060 
00001070 
00001080 
00001090 
00001100 
00001110 
00001 720 
00001130 
00001140 
00001 150 
00001160 
00001 170 
00001180 
00001 190 
00001200 
00001210 
00001220 


1 
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C IS FEN ALREADY POUND? 

00001230 

70 

If (P2N.5Q. 1 . D2) GO TO 80 

00001240 

C YES 

00001250 


WRITE (6 , DUMPS II) 

00001260 


GO TO 30 

00001270 

C 


00001280 

C CALCULATE DSDX 

00001 290 

80 

ZE1A = 2(3) *SP-Z (2} *C? 

00001300 


BETA3= BETA**3/(1. DO-BETA) 

00001310 


PZ5= 2 (2)*BETA3 

00001320 


PZ6= Z (3) *BETA3 

00001330 


DXDH= -2. D.Q*Z (2) *B3TA*CF*Z (3) *BETA*SF+PZ5*Z£TA*Z (2) 

00001340 


DXDK= Z (2)*BETA*SF-1.D0+PZ6*Z (2)*ZE?A 

00001350 


D YDH= Z ( 3J*B2TA*CF-1 . D0-PZ5*Z (3)*ZETA 

00001360 


P YDK= -2. P0*Z (3)*32TA*SF+Z (2) *BETA*CP-PZ6*Z (3)*ZETA 

00001 370 


DSDX ( 1 ) = 2.D0*Z (1 ) ** (-3) 

00001 380 


DU M 1 = 2. D0*D1*X1-D3*Y1 

00001390 


D0H2= 2.D0*D2*Yl-D3*Xl 

00001400 


DSDX (2)= D0H1*DXDH+DUM2*DYDH 

00001410 


DSDX(3) = DUM1*DXDK+DUM2*DYDK 

00001420 


D=2.D0/(1.D0*Z(4) **2*Z (5) **2) 

00001430 


DXSP- {-YSUf!*Z (5)-ZS0N) *D 

00001 Y40 


DXSQ= YSUN*Z(4)*D 

00001450 


DYSP= XSUN*Z (5)*D 

00001 460 


DYSQ= {-XSON*Z (4) +ZSUN) *D 

00001470 


DU!11= -2. DO* XI* (X1*XSUN*Y1*YSUN) 

00001480 


DUfl2= -2.D0* Y1 * {Y1*YS0N+X1*X5UN) 

00001 490 


DSDX (4) = DOM 1 *DXS? *DUM2*DY3P 

00001500 


DSDX (5) = D031*DX5C+D0fl2*r YS0 

00001510 

C 

WRITE (6 f PE4 ) 

00001520 


IS HA D= IS HA D * 1 

00001530 


IF (DSDF. LT. 0. DO) GO TO 100 

00001540 

C 


00001550 

C EXIT ANGLE AND DERIVATIVE 

00001560 


F EX= DU M 

00001570 


DO 90 J= 1 , 5 

00001580 

90 

DFE X ( J)= -D 3X ( J)/DSDF 

00001590 


GO TO 10 

00001600 

e 


00001610 

C ENTRY ANGLE AND DEBIVATXYE 

00001620 

100 

FEN=DUfI 

00001630 


DO 1 1 0 J = 1 , 5 

00001640 

1 10 

DFE N (J) = -DSDX (J)/DSDF 

00001650 


GO TO 10 

00001660 

t* 


00001670 

c 


00001680 

120 

IF ( (ISBAD.EQ.O) .OB. (ISHAD.EQ.2) ) P.ETORN 

00001690 


URITE (6,1020) ISRAD 

00001700 


ISRAD= 0 

00001710 


RETURN 

00007720 

c 


00001730 

130 

WRITE (6,1030) 11EE 

00001740 


STOP 

00001750 

c 


00001760 

c 


00001770 

1010 

FORMAT (33HO DSDF=0, ORBIT TANGENT TO SHADOW) 

00001780 

1020 

FORMAT (15H0 ERRO R--ISH AD=, 1 4) 

00001790 

1030 

FORMAT ( 4 9 HO DQRTIC HAS RETURNED WITH NUMBER OF REAL ROOTS = ,16) 

00001800 


END 

00001810 
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Subroutine SUN (SUNSS) 


Description : 

Calculates planet-sun unit vector in the equinoctial frame and 
the distance in A.U.'s divided by planet's seraimajor axis. 

Argument List : 

T, | 

T 
Z 

Common areas: 

SOL/ RS (3) , R 

TERRA/£§, £g, W, ENg, Ag, ECLMTX(3,3) , EOUMTX(3,3) 

Called by : 

FUNCT, OUTP 


Time 

State and costate 


i 
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C EUN/SUNSS 

c 

c 

C THIS PROGRAM CALCULATES THE PLANET TO SUN DIRECTION AND 
C DISTANCE FOR A GIVEN TIME, OUTPUT IN THE EQUINOCTIAL 
C COORDINATE FRAME* USED WITH PLANETS. 

C INPUT 

C Z--10 VECTOR C? EQ. OK AND COST AT E ( NOT USED) 

C AE- -0 KBIT SESIMAJOR AXIS 

C EC"-CR3IT ECCENTRICITY 

2 W--ARG. OF PERIH. 

C EN2 — MEAN ORBITAL NOTION 

C AN — MEAN ANOMALY AT BEGINNING OF TRAJECTORY (TO) 

C ECLMTX-- MATRIX FCR CONVERSION TO ECLIPTIC FRAME 

C EQUMIX— MATRIX ?0R CONVERSION TO EQUATORIAL PRAHE 

C T--PRSS2NT TIME 

C OUTPUT 

C RS — UNIT VECTOR FROM PLANET TO SUN, EQUINOCTIAL COORD. 

C R — DISTANCE FROM PLANET TO SUN AT TIME T 

C 
C 

c 

SUBROUTINE SUN (T, Z) 

C 

IMPLICIT REAL*8 (A-H,0-Z) 

C 

COMMON /SOL / RS (3) , E 

COMMON/TERRA/ AE , EC, W , EN2, AN , ECLKT2 (3,3) , E Q □ M TX (3,3) 

C 

- DIMENSION RSI (3) , CM (3,3) ,Z (10) ,RS2 (3) ,TMP (3) 


MEAN ANOMALY AT TIME T 

AA= AN+3NE*T 

TRUE ANOMALY --CORRECT THRU ECCENTRICITY CUBED 

F=AA+ (2- D0*EC-. 25D0 *EC ** 3) *D S IN (A A) + 1 .25D0*SC**2*DSIH (2. DO* A A) 
1 +1.083333333333D0*3C**3*D5IN(3.D0*AA) 

E=F+8 


DISTANCE BETWEEN PLANET AND SUN IN PLANET RADII 
R= p.DO-EC**2)/(1.DO + EC*DCOS (P) ) 

CALCULATE UNIT VECTOR TO SUN, EARTH ECLIPTIC FRAME 
IMP ( 1 ) =-DC0S (B) 

IMP ( 2 ) =~DSIN (B) 

TMP (3) =0. DO 

DO 3 1=1,3 
RS2 < ■ 0> DO 
DO 3 0=1,3 

RS2 ( I )= R 3 2 (I>ECLMTZ ( J , I) *TMP (J) 

CALCULATE UNIT VECTOR TO SUN, EQUATOEIAL FRAME 
DO 5 1=1,3 
RSI (I>O.D0 
DO 5 0=1,3 

RSI (I>RS1 (I>EQUKTX ( 0,1) *RS2 (J) 


TRANSFORM TO EQUINOCTIAL COORD. 


009000 10 
00900020 
00000030 
00000040 
00000050 
00000060 
00000070 
00000080 
00000090 
00090100 
00000110 
00000 120 
00000 130 
00000 149 
OOOOO 150 
00000 160 
00000 1 70 
00000 180 
00000190 
00000200 
00000210 
00000220 
00000230 
00000240 
00000 250 
00000260 
00000270 
00000280 
00000290- 
00000300 
00000310 
00000320 
00000330 
00000340 
00000350 
00000360 
00000370 
00000380 
00000390 
00000400 
00000410 
00000420 
00000430 
00000 440 
00000450 
00000460 
Or) 000470 
00000490 
00000490 
00000500 
00000510 
00000520 
00000530 
00000 540 
00000550 
00000560 
00000570 
00000580 
09000590 

oooooeno 

00000610 


90 


i 



AB= 1. DO«-Z (4} **2+Z (5) **2 

CM (1,1) = (l-DO-2 (4)**2+Z (5)**2) /AB 

CH(2,1) = 2. DO*Z (4) *Z (5) /AB 

CM (3 , 1) = -2,D0*Z (4)/AB 

CM (1 , 2} = CM (2,1) 

CM (2,2)- (1. DO+-Z (4)** 2-Z (5)**2) /AB 
CM (3,2)= 2-D0*Z(5)/AB 
CM C 1 # 3) = ~CM(3,1) 

CM (2,3)= -CM (3,2) 

CM (3,3) = (1 .DO-Z (4) **2-Z (5) **2)/AB 
DO 10 1=1,3 
RS(I)= O.DO 
DO 10 J=1 , 3 

10 RS (I)= RS (I)*ca (J,I) *RS1 (J) 

c 

RETURN 

END 


00000620 
0000063 0 
00000640 
00000 650 
00000660 
00000670 
00000680 
00000690 
00000700 
00000710 
00000720 
00000730 
00000740 
000007 50 
00000760 
00000770 
00000780 
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Subroutine TRAJ (TRAJSS) 


Description : 

The Runge-Kutta integrator is called to calculate the low thrust 
trajectory. The errors in the final conditions are calculated to send 
back to the iterator. 

Common Areas: 

XMMM/ ZLO ( 5) , STEP ( 6 ) , ZERF(6) 

TRA/TFMAX, gXQ, HgB, EW(IO) 

Z/ Z(10) , PBRZ (10) 

INT/ IPR. TDTM . TDIM2 . NIMAX 
T/TF, TO 

ELEM/ ZPO (5) , ZPF ( 5 ) 

DY/ DYDT ( 6) 

TC/NOP 

A/A, AMU, PX 

NU/XMIJ 

Called by : 

CONTL, ITER 

Calls Subroutines : 

RK4, FUNCT 
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C TRAJ/TFAJSS 


C 5CLAR SAIL 

C INCLUDES PENALTY FUNCTION. 

C 'IBIS ROUTINE SETS UP THE INPUT TO THE INTEGRATOR WHICH 
C EXTRAPOLATES THE TRAJECTQRY FROM INITIAL TIMS TO 
C FINAL TIME. IT ALSO EVALUATES THE CHANGE IN TF AND 
C THE ERROR IN THE FINAL CONDITIONS. 

C THIS PROGRAM IS CALLED BY ITER OR BY COHTROL 
C IT CALL THE SUBPROGRAM RK 4 (RONGA-KUTTA) 

C MIN J, MAX H. 

C 6 DIN. ZERF. T.C. OPTIOHS. 

C NCP= 1 --ALL 5 FINAL O.E. FIXED, -2— A, 2,1 ONLY FIXED. 

2 =3— SEMIMAJOR AXIS (C3) SPECIFIED , ESCAPE TRAJ 

C 

SUBROUTINE TRAJ 

C 

IKPLZCIT REA L*8 ( A-H , 0 -Z) , INTEGER (I-N) 

Q 

COMMON /XMHK/ZLO ( 5 ) , STEP (6) * ZERF(6) 

COMMON /TBA/TFMAX , DTO , OEB, ER(10) 

COMMON /Z/Z (10) , DERZ (10) 

COMMON /INT/IPR, TDIM , IDIM2, NIKAX 

CCHMON /T/TF, TO 

COMMON /ELEM/ZPO ( 5 ) , ZPF ( 5 ) 

COHBON /DY/DTDT (6) 

COMMON /TC/NOP 
COMMON /A/A, AMU, P I 
COMMON /NO/XNU 
C 

EXTERNAL FUNCT, OUTP 

DIMENSION PEMT (5) , AUX (1 , 10} , DSBZ1 (10) 

C . 

c 

9 IF (TF.GT.TFMAX) TF = TF MAX 
PRMT(1) = TO 

PRMT (2) = TF 
PRMT (3)= DTO 
PRET (4) = UEB 
C 

c 

C Z IS A VECTOR OF STATE AND COSTATE 

DO 10 1=1 , IDIM2 
Z [I) =ZP0 (I) 

10 Z {I + IDI M2) = ZLO (I) 

C 

C 28 ARE ERROR WEIGHTS- -INPUT TO THE INTEGRATOR 
C 

15 DO 20 1=1, IDIH 
20 DERZ (I)-%% (I) 

C 

C CALL THE R—K INTEGRATOR 

C 

CALL RK4 (PRMT, Z, DERZ, IDIH,IHLE, FUNCT, OUTP, AOX) 

IF (IHLF.GT.10) GC TO 100 
C 

C Z IS NOX THE FINAL O.E. , 

C ZERP THE ERROR IN THE FINAL CONDITIONS 

C 

H=0. DO 


00000010 
00000 020 
00000030 
00000040 
00003050 
00000060 
00000070 
00090080 
00000090 
00000100 
00000110 
00000120 
00000130 
00000 140 
00000150 
00000160 
00000170 
00000 180 
00000 190 
00000200 
00000210 
00000 220 
00000230 
OOOOO 240 
00000250 
00000260 
00000270 
00000280 
00000290 
00000300 
00000310 
00000320 
OOOOO 330 
OOOOO 340 
OOOOO 350 
OOOOO 360 
000 00370 
00000380 
00000390 
00000400 
0000041 0 
OOOOO 420 
00000430 
00000440 
00000450 
00000460 
00000 470 
00000480 
00 0 004 90 
00000500 
00000510 
00000520 
00000530 
00000540 
00000550 
00000560 
00000570 
00000580 
00000590 
00000600 
00000610 
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DO 30 1=1 , IDIH2 
30 H= H + Z (1+5) *DERZ (I) 

TF1=TF* (STEP (6)-i-1.D0) 

CALL FUNCT(TF1 ,Z,DERZ1) 

H1=0. DO 

DO 35 1=1, IDIM2 
35 H1=H1*Z (1+5) *DERZ1 (I) 

DYDT (6) = (H1-H)/ (TF1-TF) 

C PENALTY FUYCTION ETPECT * ** 

R = H-XNU/ (2. D0*Z (1)*Z (1 ) ( 1 . DO-DSQRT(Z (2)*Z (2)+Z (3)*Z (3))) 2) 

C 

ZERF (6) = H -1.00 

C 

C FINAL COYDITION OPTION BRANCH 
C 

GO TO (40,50,80) ,NOP 

c 

40 DO 45 1=1,5 

ZERF (I) = Z (I) -ZPF (I) 

45 DYDT (I) = DERZ (I) 

RETURN 


50 


60 


ZERF (4) = (Z (3)*Z (7)-Z (2)*Z(8) )*1.D-3 
ZERF (5) = (Z (5) *Z (9) -Z (4) *Z (10) ) *1. D-3 

DYDT (4) = DERZ (3) *Z (7) +Z (3) *DERZ (7) -DERZ (2) *Z (8) -Z (2) *DERZ (8) 
DYDT (4) = DYDT {4} *1. D-3 

DYDT (5) = DERZ (5)*Z (9) +Z (5)*DERZ (9) -DERZ (4) *Z ( 10) -Z (4) *DERZ (10) 
DYDT ( 5) = DYDT (5)* 1. D-3 

ZERF (1) = Z(l) - ZPF ( 1) 

DUM1= DSQRT (Z (2) **2 + Z(3)**2) 

ZERF (2) = D0H1 " ZPF (2) 

DUM2= DSQRT(Z (4) **2 + Z(5)**2) 

ZERF {3)= DDM2 " ZPF (3) 

DYDT ( 1)= DERZ (1) 

DYDT (2)= 0.D0 
DYDT (3) = 0. DO 

IF (DUM1 .GT. 1 .D-1 2) DYDT{2)= (Z (2) *EERZ (2) + Z ( 3) *DERZ (3) ) /DUM 1 
IF (D0H2.GT. 1.D-12) DYDT (3) = (Z ( 4) *D ERZ ( 4) + Z (5) *qeRZ (5)) /D0M2 


C 

C SPECIAL CASE, E=0 AH D/OH 1=0 

C 

IF (ZPF (2) .NE.O.DO) GO TO 70 

ZERP (2) = Z (2) 

ZERF (4) = Z (3) 

DYDT (2) = DERZ [2) 

DPDT (4) = DERZ (3) 

70 IF (ZPF (3) .NE.O.DO) RETURN 
ZERP (3)= Z {4} 

ZEHF (5)= Z (5) 

DYDT (3) = DERZ (4) 

DYDT (5)= DERZ (5) 

RETURN 


C 

80 ZERF (1)= Z (1) -ZPF (1) 

DYDT (1) = DERZ (1) 

DO 85 1=2,5 
ZERF (I)= Z { I + 5 ) *1 . D-3 
85 DYDT (I)= DERZ (I+5)*1. D-3 
RETURN 
C 


00000620 
00030630 
00000640 
00000650 
00000660 
00000670 
00000680 
00000690 
00000 700 
00000710 
00000720 
00000730 
00000740 
00000750 
00000760 

0 0000 7 70 
00000780 
00000790 
00000 800 
00000R10 
00000820 
00000830 
00000840 
00000850 
00000 860 
00000870 
00000880 
00000890 
00000900 
00000910 
00000920 
00000930 
00000940 
0000095 0 
00000960 
00000970 
00000980 
00000990 
00001 000 
00001010 
00001020 
00001030 
00001040 
00001 050 
00001060 
00001070 
00007080 
00001090 
00001 100 

00001 1 10 
00001 120 
00001 130 
00001 140 
00001 150 
00001 160 
00001 1 70 
00001 1 80 
00001 190 
00001200 
00001210 
00001220 





j 



100 IF (XHLF-EQ.il) WRITS (6,1000) 
IF (IHLF.EQ. 12) WRITE (5,1001) 
If (IHLF.EQ- 1 3) WRITE (6,1002) 

STOP 


00001 230 
00001 240 
00001250 
00001260 
00001 270 


1000 FORMAT (68H0 THE NUM3ZE OF BISECTIONS OF THE ORIGINAL INCREMENT HA00001 280 

IS EXCEEDED 10) 00031290 

1001 FORMAT (27H0 INITIAL INCREMENT IS ZERO) 00001300 

1002 FORM A' 1 ’ ( 5 4 HO INITIAL INCREMENT HAS WRONG SIGN OR BOUNDS ARE WRONG) 00001 310 

pwn * 00001320 
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Subroutine VCP 


Description: 

Calculates a vector cross product, W = U x V 

Argument List: 

12 First 3-vector 

V Second 3- vector 

K. Cross product result, U x V 
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C VCP 

THIS SUBROUTINE CALCULATES THE VECTOR CROSS PRODUCT 0 X V 
AND FETURNS PRODUCT AS V 

SUBROUTINE VCP(U, V,W) 

DOUBLE PRECISION U,V,W 
DIMENSION U (3) , V (3) ,V (3) 

R ( 1 )= -0(3)*V (2)*XJ (2)*V( 3) 

K (2)= 0(3)*V(1)-0(1)*V(3; 

W (3) = -U (2)*V(1)+U (1)*V (2) 

EE TURN 

END 


00000010 
000QU0 20 
OOOOO 030 
00000040 
00000050 
09000050 
00000070 
OOOOO 080 
00000090 
00300 100 
05000110 
0C000 120 
00000130 
00000 140 
00000150 
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